format: text

text:
# Dynamics and Rheology of a Supercooled Polymer Melt in Shear Flow

## I Introduction

The dynamics and rheology of glassy polymers are known to be very complicated and are still not well understood. We summarize salient features of such systems below.

First, in the linear response regime, thermal relaxations of the chain conformations occur from microscopic to macroscopic time scales, as revealed in measurements of stress and dielectric responses. In a relatively early stage, the stress relaxation function $`G(t)`$, which describes linear response to small shear deformations, can be fitted to the Kohlrausch-Williams-Watts (KWW) form,

$$G_\mathrm{G}(t)=G_0\mathrm{exp}[(t/\tau _\mathrm{s})^c],$$
(1.1)
after a microscopic transient time $`t_{\mathrm{tra}}`$. The time $`\tau _\mathrm{s}(\tau _{\mathrm{tra}})`$ is of the order of the structural $`\alpha `$ relaxation time $`\tau _\alpha `$ (to be defined in (3.8) below), which grows dramatically as the temperature $`T`$ is lowered towards the glass transition temperature $`T_\mathrm{g}`$. When the time $`t`$ considerably exceeds $`\tau _\mathrm{s}`$, the relaxation of the chain conformations is relevant and is well described by the Rouse or reptation dynamics, depending on $`N<N_\mathrm{e}`$ or $`N>N_\mathrm{e}`$, respectively. Here $`N`$ is the polymerization index and $`N_\mathrm{e}`$ is that between entanglements on a chain. For short chain systems with $`N<N_\mathrm{e}`$, the overall behavior in the time region $`tt_{\mathrm{tra}}`$ may be expressed as

$$G(t)=G_\mathrm{G}(t)+G_\mathrm{R}(t).$$
(1.2)
The $`G_\mathrm{R}(t)`$ is the stress time-correlation function in the Rouse model whose terminal relaxation time $`\tau _\mathrm{R}`$ is of order $`N^2\tau _\alpha `$. For entangled chain systems with $`N>N_\mathrm{e}`$, the KWW function in (1.1) is followed by the power-law decay,

$$G(t)e^1G_0(t/\tau _\mathrm{s})^\nu $$
(1.3)
with $`\nu 0.5`$ until the rubbery plateau $`G(t)G_N^{(0)}`$ is reached, where $`G_N^{(0)}`$ assumes the modulus $`nk_\mathrm{B}T/N_\mathrm{e}`$ of entangled polymers with $`n`$ being the bead number density. Ultimately, $`G(t)`$ follows the reptation relaxation $`G_{\mathrm{rep}}(t)`$ on the time scale of a very long reptation time $`\tau _{\mathrm{rep}}`$. These hierarchical relaxations arise from rearrangements of jammed atomic configurations and subsequent evolution of chain conformations. They also give rise to the corresponding characteristic behaviors in the frequency-dependent shear modulus $`G^{}(\omega )i\omega _0^{\mathrm{}}𝑑te^{i\omega t}G(t)`$, depending on the frequency $`\omega `$ relative to the inverse characteristic times introduced so far.

Second, in the nonlinear response regime, glassy fluids generally exhibit highly viscous non-Newtonian flow close to (but above) $`T_\mathrm{g}`$ even if they are low-molecular-weight fluids. In such fluids, if $`\dot{\gamma }>\tau _\alpha ^1`$, atomic rearrangements are induced not by thermal agitations but by externally applied shear. In chain systems without entanglements, on one hand, shear thinning occurs at sufficiently high (but sometimes unrealistically large) shear rates due to chain elongation. In entangled polymers, on the other hand, shear thinning occurs at very small shear larger than $`\tau _{\mathrm{rep}}^1`$, where disentanglements are induced by shear. Thus supercooled chain systems are most easily driven into a nonlinear response regime even by extremely small shear, though the crossover shear stress from linear to nonlinear regimes may not be very small. Furthermore, in glassy fluids below $`T_\mathrm{g}`$, plastic deformations are often induced in the form of large-scale shear bands above a yield stress (corresponding to a few $`\%`$ strain).  It is of great importance to understand how these nonlinear effects occur depending on $`\dot{\gamma }`$, $`T`$, and $`N`$.

Third, in rheological experiments on polymers, use has been made of the stress-optical relation between the deviatoric (anisotropic) parts of the dielectric tensor $`ϵ_{\alpha \beta }`$ (at optical frequencies) and the average stress tensor $`\sigma _{\alpha \beta }`$. In shear flow with mean velocity $`\dot{\gamma }y`$ in the $`x`$ direction, it is expressed as

$$ϵ_{xy}=C_0\sigma _{xy},ϵ_{\alpha \alpha }ϵ_{\beta \beta }=C_0(\sigma _{\alpha \alpha }\sigma _{\beta \beta }),$$
(1.4)
where $`C_0`$ is called the stress-optical coefficient. For melts this relation excellently holds at relatively high $`T`$ ($`>T_\mathrm{g}`$) for general time-dependent nonlinear shear deformations. If measurements are made in steady states, it holds even close to $`T_\mathrm{g}`$. For its validity, we need to require that the form contribution to $`ϵ_{\alpha \beta }`$ is negligible as compared to the intrinsic contribution and that the glassy part of the stress is negligible as compared to the usual entropic part. Thus it is violated when the form part is relevant such as in polymer solutions close to the demixing critical point or when measurements are made in transient states close to $`T_\mathrm{g}`$. In the latter case, as is evident from the enhancement of $`\tau _\mathrm{s}`$ in (1.1), the glassy part of the stress is dominant for relatively rapid deformations.

While the predictive power of analytic theories in polymer science is still poor, computer simulations can provide us a useful tool to investigate the microscopic origins of experimentally observed macroscopic phenomena. In quiescent states, diffusive motions in supercooled melts have been extensively studied using molecular dynamics (MD) and Monte Carlo simulations. In another application, nonequilibrium molecular dynamics (NEMD) simulations have been useful to investigate chain deformations and rheology in flow. In particular, Kröger et al. studied the molecular mechanisms of the violations of stress-optical behavior for a melt consisting $`M=260`$ chains with bead number $`N=30`$ after application of elongational flow.

In this paper, we will present results of very long MD simulations to study linear and nonlinear dynamics of a supercooled polymer melt in the absence and presence of shear flow. Long simulation times are needed to calculate the terminal relaxation of $`G(t)`$, which has not yet been undertaken in the literature. As a new finding, we will show that each chain in our melt system is changing its orientation (tumbling) randomly in shear flow. Use will be made of techniques and concepts introduced in our previous papers on supercooled binary mixtures under shear flow. Some of our results were already published elsewhere.

## II Model and simulation method

Our system is composed of $`M=100`$ chains with $`N=10`$ beads confined in a cubic box with length $`L=10\sigma `$ and volume $`V=L^3=10^3\sigma ^3`$. The number density is fixed at $`n=NM/V=1/\sigma ^3`$, which results in severely jammed configurations at low $`T`$. All the bead particles interact with a truncated Lennard-Jones potential defined by

$$U_{\mathrm{LJ}}(r)=4ϵ\left[\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6\right]+ϵ(r<2^{1/6}\sigma ).$$
(2.1)
The right hand side is minimum at $`r=2^{1/6}\sigma `$ and the potential is truncated for larger $`r`$ ($`U_{\mathrm{LJ}}(r)=0`$ for $`r>2^{1/6}\sigma `$). By using the repulsive part of the Lennard-Jones potential only, we may prevent spatial overlap of the particles. Consecutive beads on each chain are connected by an anharmonic spring of the form,

$$U_\mathrm{F}(r)=\frac{1}{2}k_cR_0^2\mathrm{ln}[1(r/R_0)^2]$$
(2.2)
with $`k_c=30ϵ/\sigma ^2`$ and $`R_0=1.5\sigma `$. In our simulation the bond lengths $`b_j^k|𝑹_j^k𝑹_{j+1}^k|`$ ($`1jN1)`$ between consecutive beads on the same chain $`k`$ were very close to the minimum distance $`b_{\mathrm{min}}0.96\sigma `$ of the sum $`U_{\mathrm{LJ}}(r)+U_\mathrm{F}(r)`$. The deviations $`b_j^kb_{\mathrm{min}}`$ were only on the order of a few $`\%`$ of $`b_{\mathrm{min}}`$ for any $`T`$ and $`\dot{\gamma }`$ realized in our study.

Microscopic expressions for physical quantities can be expressed in terms of the momentum and position vectors of the $`j`$-th bead on the $`k`$-th chain, $`𝑹_j^k`$ and $`𝒑_j^k`$, where $`j=1,\mathrm{},N`$ and $`k=1,\mathrm{},M`$. For example, the space integral of the microscopic stress tensor reads

$`\mathrm{\Pi }_{\alpha \beta }^\mathrm{T}(t)`$ $`=`$ $`{\displaystyle \frac{1}{m}}{\displaystyle \underset{k=1}{\overset{M}{}}}{\displaystyle \underset{j=1}{\overset{N}{}}}p_{j\alpha }^kp_{j\beta }^k{\displaystyle \underset{\mathrm{all}\mathrm{pairs}}{}}U_{\mathrm{LJ}}^{}(\xi ){\displaystyle \frac{\xi _\alpha \xi _\beta }{\xi }}`$ (2.4)
$`{\displaystyle \underset{k=1}{\overset{M}{}}}{\displaystyle \underset{j=1}{\overset{N1}{}}}U_\mathrm{F}^{}(\xi ){\displaystyle \frac{\xi _\alpha \xi _\beta }{\xi }},`$
where $`m`$ is the mass of a bead, $`U_{\mathrm{LJ}}^{}(\xi )=dU_{\mathrm{LJ}}(\xi )/d\xi `$, and $`U_\mathrm{F}^{}(\xi )=dU_\mathrm{F}(\xi )/d\xi `$. Here $`𝝃=(\xi _x,\xi _y,\xi _z)`$ in the right hand side represents the relative vector $`𝑹_j^k𝑹_j^{}^k^{}`$ between the two beads, $`𝑹_j^k`$ and $`𝑹_j^{}^k^{}`$, in the second term and the relative vector $`𝑹_j^k𝑹_{j+1}^k`$ between the two consecutive beads, $`𝑹_j^k`$ and $`𝑹_{j+1}^k`$, of the same chain in the third term. To avoid cumbersome notation, we will write the bead positions simply as $`𝑹_j`$ $`(j=1,\mathrm{},N)`$ suppressing the index $`k`$. When they will appear in the statistical averages $`\mathrm{}`$, the average over all the chains $`_{k=1}^M(\mathrm{})/M`$ will be implied even if not written explicitly. Furthermore, it is convenient here to introduce the usual notation $`\sigma _{\alpha \beta }`$ for the stress tensor by

$$\frac{1}{V}\mathrm{\Pi }_{\alpha \beta }^\mathrm{T}=p\delta _{\alpha \beta }\sigma _{\alpha \beta },$$
(2.5)
where $`p`$ is the pressure and the second term is deviatoric. The $`\sigma _{\alpha \beta }`$ has already appeared in the stress-optical relation (1.4).

We will hereafter measure space and time in units of $`\sigma `$ and $`\tau _0(m\sigma ^2/ϵ)^{1/2}`$. The temperature $`T`$ will be measured in units of $`ϵ/k_\mathrm{B}`$. The original units will also be used when confusion may occur. Our simulations cover normal ($`T=1.0`$) and supercooled ($`T=0.2`$) states with and without shear flow ($`\dot{\gamma }=0`$, $`10^4`$, $`10^3`$, $`10^2`$, and $`10^1`$). Simulation data were taken after very long equilibration periods ($`10^2\tau _R5\times 10^6`$ at $`T=0.2`$) so that no appreciable aging (slow equilibration) effect was detected in the course of taking data in various quantities such as the pressure or the density time-correlation functions. (i) In quiescent cases, we impose the micro-canonical condition and integrate the Newton’s equations of motion,

$$\frac{d}{dt}𝑹_j=\frac{1}{m}𝒑_j,\frac{d}{dt}𝒑_j=𝒇_j,$$
(2.6)
where $`𝒇_j`$ is the force acting on the particle $`j`$ due to the potentials. Integration was performed with time increment $`\mathrm{\Delta }t=0.005`$ under the periodic boundary condition. Long time simulations of order $`10^2\tau _\mathrm{R}`$, which corresponds to $`10^9`$ MD steps for $`T=0.2`$, were performed. In the previous simulations, however, the integrated times did not much exceed $`\tau _\mathrm{R}`$ in supercooled states. (ii) In the presence of shear, rewriting the momentum deviations $`𝒑_jm\dot{\gamma }Y_j𝒆_x`$ from the mean flow as $`𝒑_j`$, we integrated the so-called SLLOD equations of motion,

$`{\displaystyle \frac{d}{dt}}𝑹_j`$ $`=`$ $`{\displaystyle \frac{1}{m}}𝒑_j+\dot{\gamma }Y_j𝒆_x,`$ (2.7)
$`{\displaystyle \frac{d}{dt}}𝒑_j`$ $`=`$ $`𝒇_j\dot{\gamma }p_{y_j}𝒆_x\widehat{\zeta }𝒑_j,`$ (2.8)
where $`𝒆_x`$ is the unit vectors in the $`x`$ (flow) direction, $`𝑹_j=(X_j,Y_j,Z_j)`$, and $`𝒑_j=(p_{xj},p_{yj},p_{zj})`$. The friction coefficient $`\widehat{\zeta }`$ was set equal to

$$\widehat{\zeta }=\underset{j}{}(𝒇_j𝒑_j\dot{\gamma }p_{xj}p_{yj})/\underset{j}{}𝒑_j^2.$$
(2.9)
The temperature $`T`$ $`(_j𝒑_j^2/mNM`$) could then be kept at a desired value. The time increment was $`\mathrm{\Delta }t=0.0025`$. After an equilibration run in a quiescent state for $`t<0`$, we gave all the particles the average flow velocity $`\dot{\gamma }Y_j𝒆_x`$ at $`t=0`$ and then imposed the Lee-Edwards boundary condition to maintain the shear flow. Steady sheared states were realized after transient relaxations.

## III Dynamics in quiescent states

Although it is highly nontrivial, it has been confirmed by computer simulations that the single-chain near-equilibrium dynamics in unentangled melts can be reasonably well described by (or mapped onto) the simple Rouse model. In the Rouse dynamics, the relaxation time of the $`p`$-th mode of a chain is expressed in terms of a friction coefficient $`\zeta `$ and a segment length $`b`$ as

$$\tau _p=\zeta b^2/[12k_\mathrm{B}T\mathrm{sin}^2(\pi p/2N)],$$
(3.1)
where $`1pN1`$. The Rouse relaxation time $`\tau _\mathrm{R}`$ is the slowest relaxation time,

$$\tau _\mathrm{R}\tau _1N^2\zeta b^2/(3\pi ^2k_\mathrm{B}T).$$
(3.2)
The segment length $`b`$ in the corresponding Rouse model may be related to the variance of the end-to-end vector of a chain $`𝑷𝑹_N𝑹_1`$ in our microscopic model by

$$|𝑷|^2=b^2(N1).$$
(3.3)
As a result, $`b`$ is dependent on $`T`$ but its dependence turns out to be weak as $`b=1.17`$, $`1.18,`$ $`1.19`$ for $`T=1.0`$, $`0.4`$, $`0.2`$, respectively. Note that $`b`$ is larger than the minimum distance $`b_{\mathrm{min}}0.96`$ of the bond-potential. Let us consider the time-correlation function of $`𝑷(t)`$,

$$C(t)=𝑷(t+t_0)𝑷(t_0)/|𝑷|^2,$$
(3.4)
which is normalized such that $`C(0)=1`$. Here $`C(t)`$ should be independent of the initial time $`t_0`$ in steady states in the limit of large system size. However, our system is not very large, so we took the average over the initial time $`t_0`$. This statistical averaging will not be mentioned hereafter in showing our MD results of time-correlation functions. In the Rouse dynamics $`C(t)`$ is calculated as

$$C_\mathrm{R}(t)=\frac{2}{(N1)N}\underset{\mathrm{odd}p}{}\mathrm{cot}^2\left(\frac{\pi p}{2N}\right)e^{t/\tau _p},$$
(3.5)
where the summation is over odd $`p`$ but the first term ($`p=1`$) is dominant in the whole time region (so we may determine $`\tau _\mathrm{R}`$ by $`C(\tau _\mathrm{R})=e^1`$). As shown in Fig.1, our MD data of $`C(t)`$ can be excellently fitted to $`C_\mathrm{R}(t)`$. The $`\tau _\mathrm{R}`$ thus determined increases drastically with lowering $`T`$ as $`\tau _\mathrm{R}=250`$, $`1800`$, and $`6\times 10^4`$ for $`T=1.0`$, $`0.4`$, and $`0.2`$, respectively. In the previous simulations on nonentangled polymer melts,  numerical results were consistent with the Rouse dynamics for small $`p`$ (large-scale motions), but deviations are enhanced for large $`p`$ (small-scale motions) in supercooled states. Furthermore, we here give the expression for the stress relaxation function in the Rouse model,

$$G_\mathrm{R}(t)=\frac{nk_\mathrm{B}T}{N}\underset{p=1}{\overset{N1}{}}\mathrm{exp}(2t/\tau _p),$$
(3.6)
which is equal to $`nk_\mathrm{B}T(N1)/N`$ at $`t=0`$ and decays as $`nk_\mathrm{B}TN^1\mathrm{exp}(2t/\tau _\mathrm{R})`$ for $`t\tau _\mathrm{R}`$. Since $`G(t)`$ is much larger than $`G_\mathrm{R}(t)`$ in the relatively short time region $`t<\tau _s`$, they can coincide only in the late stage.

In Fig.2, we show the van Hove self-correlation function,

$$F_q(t)=\frac{1}{N}\underset{j=1}{\overset{N}{}}\mathrm{exp}[i𝒒\mathrm{\Delta }𝑹_j(t)],$$
(3.7)
where $`q=2\pi `$, $`\mathrm{\Delta }𝑹_j(t)=𝑹_j(t+t_0)𝑹_j(t_0)`$ is the displacement vector of the $`j`$-th bead in the time interval $`[t_0,t_0+t]`$. The peak wave number of the static structure factor is given by $`q2\pi `$. We define the $`\alpha `$ relaxation time $`\tau _\alpha `$ from the condition,

$$F_q(\tau _\alpha )=e^1(q=2\pi ).$$
(3.8)
As has been reported in the literature, $`\tau _\alpha `$ increases drastically with lowering $`T`$. In our case, we obtained $`\tau _\alpha =0.91`$, $`5.8`$, and $`310`$ for $`T=1.0`$, $`0.4`$, and $`0.2`$, respectively. At $`T=0.2`$, where the particle motions are considerably jammed, $`F_q(t)`$ exhibits a two-step relaxation and may be excellently fitted to the stretched exponential decay ($`\mathrm{exp}[(t/\tau _\alpha )^{0.64}]`$) for $`t10`$. Thus our system at $`T=0.2`$ has characteristic features of a supercooled state, although its melting temperature is unknown. We find $`\tau _\alpha 10^2\zeta b^2/k_\mathrm{B}T`$ and $`\tau _\mathrm{R}N^2\tau _\alpha `$ at any $`T`$. Particularly, for $`T=0.2`$, we obtain

$$\tau _\alpha 0.017\zeta b^2/k_\mathrm{B}T,\tau _\mathrm{R}1.9N^2\tau _\alpha .$$
(3.9)
The friction coefficient $`\zeta `$ in the mapped Rouse model grows strongly as $`T`$ is lowered in supercooled states.

Now we discuss the linear viscoelastic behavior in supercooled states. In terms of $`\mathrm{\Pi }_{xy}^\mathrm{T}(t)`$ in (2.3), the stress relaxation function $`G(t)`$ is written as

$$G(t)=\mathrm{\Pi }_{xy}^\mathrm{T}(t+t_0)\mathrm{\Pi }_{xy}^\mathrm{T}(t_0)/k_\mathrm{B}TV.$$
(3.10)
In Fig.3, we show numerical data of $`G(t)`$, where the average over the initial time $`t_0`$ was taken but the data become noisy at very large $`t\tau _\mathrm{R}`$. In the very early stage $`t1`$, $`G(t)`$ oscillates rapidly due to the vibrations of the bond vectors $`𝒃_j=𝑹_j𝑹_{j+1}`$. The initial value $`G(0)`$ takes a large value ($`100`$ in units of $`ϵ/\sigma ^3`$) nearly independent of $`T`$. For $`T=0.2`$, $`G(t)`$ can be nicely fitted to the the stretched exponential form (1.1) with $`G_05`$, $`\tau _\mathrm{s}=900.33\tau _\alpha `$, and $`c=0.5`$ in the time region $`1t10\tau _\mathrm{s}`$. For $`t50\tau _\mathrm{s}`$ it approaches the Rouse stress relaxation function $`G_\mathrm{R}(t)`$ in (3.5). The zero-frequency Newtonian viscosity is given by $`\eta (0)=_0^{\mathrm{}}𝑑tG(t)`$, so it consists of the glassy (monomeric) part,

$$\eta _\mathrm{G}=_0^{\mathrm{}}𝑑tG_\mathrm{G}(t)10\tau _\mathrm{s},$$
(3.11)
and the Rouse (polymeric) part,

$$\eta _\mathrm{R}=_0^{\mathrm{}}𝑑tG_\mathrm{R}(t)0.808TN^1\tau _\mathrm{R}.$$
(3.12)
The ratio $`\eta _\mathrm{G}/\eta _\mathrm{R}`$ is of order $`1/TN`$. They are of the same order in the present case of $`N=10`$ and $`T=0.2`$. However, we should have $`\eta _\mathrm{G}\eta _\mathrm{R}`$ for much larger $`N`$.

To examine the orientation of the bonds, we consider the orientational tensor,

$$Q_{\alpha \beta }(t)=\frac{1}{M}\underset{\mathrm{chain}}{}\frac{1}{N1}\underset{j=1}{\overset{N1}{}}\frac{b_{j\alpha }}{b_0}\frac{b_{j\beta }}{b_0},$$
(3.13)
where $`b_0^1𝒃_j`$ are the normalized bond vectors since $`|𝒃_j|b_0`$ as stated below (2.2). Notice that in the Rouse model the space integral of the entropic stress tensor is given by the expression $`\sigma _{\alpha \beta }^\mathrm{b}(3k_\mathrm{B}Tb_0^2/b^2)Q_{\alpha \beta }`$, where $`b`$ determined by (3.3) appears instead of $`b_0`$. To compare our microscopic system and the simplified Rouse model, we calculated the time-correlation function,

$$G_\mathrm{b}(t)=\sigma _{xy}^\mathrm{b}(t+t_0)\sigma _{xy}^\mathrm{b}(t_0)/k_\mathrm{B}TV,$$
(3.14)
by integrating (2.5). As shown in Fig.4, $`G_\mathrm{b}(t)`$ is fairly close to $`G_\mathrm{R}(t)`$ in (3.5) from the Rouse model. In particular, for $`t0.1\tau _\mathrm{R}`$, we find $`G(t)G_\mathrm{b}(t)G_\mathrm{R}(t)`$.

## IV Steady state behavior in shear flow

In Fig.5, we display the steady-state viscosity $`\eta (\dot{\gamma })\sigma _{xy}/\dot{\gamma }`$ obtained at $`T=0.2,0.4,`$ and $`1`$, where the time average of the stress was taken. The crossover shear rate from Newtonian to shear-thinning behavior is given by $`\tau _\mathrm{R}^1N^2\tau _\alpha ^1`$. We may introduce the Weisenberg number $`Wi`$ by

$$Wi=\dot{\gamma }\tau _\mathrm{R}\dot{\gamma }N^2\tau _\alpha .$$
(4.1)
In the non-Newtonian regime, we have $`Wi>1`$. The shear stress at the crossover is of order $`nk_\mathrm{B}TN^1`$, which is the elastic modulus of the Rouse model. The horizontal arrows indicate the linear Rouse viscosity $`\eta _\mathrm{R}`$ in (3.12), while the vertical arrows indicate the points at which $`\dot{\gamma }=\tau _\mathrm{R}^1`$. In particular, the curve of $`T=0.2`$ may be fitted to

$$\eta (\dot{\gamma })\dot{\gamma }^\nu $$
(4.2)
with $`\nu 0.7`$ for $`\dot{\gamma }\tau _\mathrm{R}1`$. The $`\eta (\dot{\gamma })`$ becomes insensitive to $`T`$ for very high shear ($`Wi1`$). In MD simulations of short chain systems in normal liquid states, similar shear-thinning has been reported, where the crossover shear is much higher, however. In MD simulations of supercooled simple binary mixtures, shear-thinning becomes apparent for $`\dot{\gamma }\tau _\alpha ^1`$, where $`\tau _\alpha `$ grows at low $`T`$ and which is consistent with (4.1) if we set $`N=1`$.

To demonstrate the stress-optical law in steady states, we show steady-state data of $`\sigma _{xy}/T`$ vs $`Q_{xy}`$ for $`T=0.2,0.4,`$ and $`1`$ in Fig.6. If the electric polarization tensor of a bead is uniaxial along the bond direction, the deviatoric part of the dielectric tensor is proportional to that of the tensor $`Q_{\alpha \beta }`$ in (3.13). In accord with the experiment, our data collapse onto a universal curve independent of $`T`$ both in the linear $`(Q_{xy}0.05)`$ and nonlinear $`(Q_{xy}0.05)`$ regimes.

We next consider anisotropy of chain conformations in shear flow. In Fig.7(a), we plot the $`x`$-$`y`$ cross-section ($`z=0`$) of the steady state bead distribution function,

$$g_\mathrm{s}(𝒓)=\frac{1}{N}\underset{j=1}{\overset{N}{}}\delta (𝑹_j𝑹_\mathrm{G}𝐫),$$
(4.3)
where $`\dot{\gamma }=10^4`$, $`T=0.2`$, and $`𝑹_\mathrm{G}=N^1_{j=1}^N𝑹_j`$ is the center of mass of a chain. In Fig.7(b), we also plot the structure factor in the $`q_x`$-$`q_y`$ plane ($`q_z=0`$),

$$S(𝒒)=\frac{1}{N^2}\underset{i,j=1}{\overset{N}{}}\mathrm{exp}[i𝒒(𝑹_i𝑹_j],$$
(4.4)
which is proportional to the scattering intensity from labeled chains in shear. In these figures $`\theta `$ is the relative angle of the ellipses with respect to the $`y`$ (shear gradient) direction. These figures demonstrate high elongation of the chains for $`\dot{\gamma }>\tau _\mathrm{R}^1`$. As will be shown in Fig.8 below, they almost saturate into the shapes shown in Fig.7 once $`\dot{\gamma }`$ exceeds $`\tau _\mathrm{R}^1`$.

Let us define the tensor,

$`I_{\alpha \beta }`$ $`=`$ $`{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{N}{}}}(R_{i\alpha }R_{j\alpha })(R_{i\beta }R_{j\beta })`$ (4.5)
$`=`$ $`{\displaystyle \frac{2}{N}}{\displaystyle \underset{j=1}{\overset{N}{}}}(R_{j\alpha }R_{\mathrm{G}\alpha })(R_{j\beta }R_{\mathrm{G}\beta }).`$ (4.6)
For small $`𝒒=(q_x,q_y,0)`$, $`S(𝒒)`$ is expanded as

$`S(𝒒)`$ $`=`$ $`1{\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha ,\beta =x,y}{}}I_{\alpha \beta }q_\alpha q_\beta +\mathrm{}`$ (4.7)
$`=`$ $`1{\displaystyle \frac{1}{2}}a_1^2(𝒒𝒆_1)^2{\displaystyle \frac{1}{2}}a_2^2(𝒒𝒆_2)^2+\mathrm{},`$ (4.8)
where $`\{𝒆_1,𝒆_2\}`$ and $`\{a_1^2,a_2^2\}`$ are the unit eigenvectors and eigenvalues of the tensor $`I_{\alpha \beta }`$ ($`\alpha ,\beta x,y`$). The two lengths $`a_1`$ and $`a_2`$ correspond to the shorter and longer radii in the principal axes of the ellipses. In terms of $`\theta `$, we have $`𝒆_1=(\mathrm{sin}\theta ,\mathrm{cos}\theta )`$ and $`𝒆_2=(\mathrm{cos}\theta ,\mathrm{sin}\theta )`$ in the $`x`$-$`y`$ plane. In Fig.8, we display $`\mathrm{tan}\theta =e_{1y}/e_{1x}`$, $`1a_1/a_2`$, and the $`xy`$ component of the alignment tensor $`Q_{xy}`$ in (3.13). All these quantities represent the degree of deformations of chain conformations in shear flow. They are insensitive to $`T`$ if plotted vs $`\dot{\gamma }\tau _\mathrm{R}`$. For $`\dot{\gamma }\tau _\mathrm{R}1`$, $`\mathrm{tan}\theta `$ is close to $`1`$ ($`\theta 45^{}`$) and both $`1a_1/a_2`$ and $`Q_{xy}`$ linearly increase with increasing $`\dot{\gamma }\tau _\mathrm{R}`$. For $`\dot{\gamma }\tau _\mathrm{R}1`$, these quantities saturate into limiting values. At $`T=0.2`$, they are

$$\theta 80^{},a_1/a_20.3,Q_{xy}0.1.$$
(4.9)
These are consistent with $`Q_{xy}\mathrm{sin}\theta \mathrm{cos}\theta `$.

## V Transient viscoelastic behavior

In Fig.9, we plot the viscosity growth function $`\sigma _{xy}(t)/\dot{\gamma }`$ after application of shear $`\dot{\gamma }`$ at $`t=0`$ for $`T=0.2`$. The system was at rest for $`t<0`$. In the initial stage $`\dot{\gamma }t0.1`$, it evolves following the linear viscoelastic growth,

$$\frac{1}{\dot{\gamma }}\sigma _{xy}(t)_0^t𝑑t^{}G(t^{}).$$
(5.1)
In the nonlinear regime, $`\sigma _{xy}(t)/\dot{\gamma }`$ tends to the non-Newtonian viscosity $`\eta (\dot{\gamma })`$. As a guide, we also display the linear growth function $`_0^t𝑑t^{}G_\mathrm{R}(t^{})`$ in the Rouse model. In the very early time region $`1t\tau _\alpha `$, the growth ($`G_0t`$) is much larger than the Rouse initial growth ($`k_\mathrm{B}Tt`$). We can see a rounded peak for $`\dot{\gamma }=0.1`$ before approach to the steady state in Fig.5. More pronounced overshoot was already reported at high shear in MD simulation of much longer unentangled alkane chains (C<sub>100</sub>H<sub>202</sub>).

In experiments, the stress-optical relation is transiently violated at low $`T`$ after application of elongational flow due to the enhancement of the glassy component of the stress. For shear flow, Fig.10 displays our MD results at $`T=0.2`$ after application of shear at $`t=0`$ in a stress-optical diagram, where the solid lines are the averages over ten independent runs. As time goes on, the system traces the curve of a given $`\dot{\gamma }`$, passes across the dashed curve representing the steady-state universal relation in Fig.6, and finally comes back to the dashed curve. The deviation from the steady-state curve becomes larger with increasing $`\dot{\gamma }`$, as in the experiments of elongational flow.

## VI Time-correlation functions and tumbling in shear flow

In Fig.11, we show the end-to-end vector correlation functions $`C(t)=𝑷(t)𝑷(0)/|𝑷|^2`$ with and without shear flow for various $`\dot{\gamma }`$ at $`T=0.2`$. For $`\dot{\gamma }0`$, it rapidly decreases, negatively overshoots, and finally approaches zero with dumped oscillation superimposed. This oscillatory behavior arises from random rotation of chains in shear flow, which is well known in dilute polymer solutions but has not been reported in polymer melts. This is more evidently seen in Fig.12, where we show time development of the $`x`$ component of the end-to-end vector $`𝑷_j=𝑹_N𝑹_1`$ of one chain for $`\dot{\gamma }=10^3`$ (a), $`10^2`$ (b), and $`10^1`$ (c) at $`T=0.2`$. The corresponding Weisenberg number (4.1) is given by $`Wi=60,600,`$ and $`6000`$, respectively. In Fig.12(d), we show chain contours projected onto the $`x`$-$`y`$ plane at the points $`18`$ indicated in Fig.12(c). When the chains change their orientation, their shapes are contracted as in the case of a single chain in solution. The average period of tumbling is about $`35/\dot{\gamma }`$ in our case.

We may introduce the van Hove time-correlation function (2.14) even in shear flow if the particle displacement vector is redefined as

$`\mathrm{\Delta }𝑹_j(t)`$ $`=`$ $`𝑹_j(t+t_0)𝑹_j(t_0)`$ (6.1)
$``$ $`\dot{\gamma }{\displaystyle _0^t}𝑑t^{}Y_\mathrm{G}(t_0+t^{})𝒆_x,`$ (6.2)
where $`Y_\mathrm{G}`$ is the $`y`$-component of $`𝑹_\mathrm{G}=N^1_{j=1}^N𝑹_j`$. From the net displacement, the first term, we have subtracted the flow-induced displacement, the second term. Figure 13 shows $`F_q(t)`$ with $`q=2\pi `$ for various $`\dot{\gamma }`$ at $`T=0.2`$. Comparison of this figure with Fig.2 suggests that applying shear is analogous to raising the temperature. This tendency was already reported for the case of supercooled binary mixtures.

We introduce the shear-dependent Rouse time $`\tau _\mathrm{R}=\tau _\mathrm{R}(T,\dot{\gamma })`$ and the $`\alpha `$ relaxation time $`\tau _\alpha =\tau _\alpha (T,\dot{\gamma })`$ by

$$C(\tau _\mathrm{R})=e^1,F_q(\tau _\alpha )=e^1.$$
(6.3)
We may then examine how shear can accelerate the motions of chains and individual beads in shear flow. Figure 14 shows $`\tau _\mathrm{R}`$ and $`\tau _\alpha `$ as functions of $`\dot{\gamma }`$ at $`T=0.2`$. In our short chain system, both $`\tau _\mathrm{R}`$ and $`\tau _\alpha `$ decrease for $`\dot{\gamma }\tau _\mathrm{R}(T,0)^1\tau _\alpha (T,0)^1`$. Our data at $`T=0.2`$ are consistent with

$$\tau _\mathrm{R}(T,\dot{\gamma })^1=\tau _\mathrm{R}(T,0)^1[1+A_\mathrm{R}\dot{\gamma }],$$
(6.4)
$$\tau _\alpha (T,\dot{\gamma })^1=\tau _\alpha (T,0)^1[1+(A_\alpha \dot{\gamma })^\mu ],$$
(6.5)
where $`A_\mathrm{R}10^4\tau _R(T,0)`$, $`A_\alpha 600020\tau _\alpha (T,0)`$, and $`\mu 0.77`$. The average tumbling period in Fig.12 is about $`4\tau _\mathrm{R}(T,\dot{\gamma })`$. For simple supercooled liquids we already introduced the van Hove time-correlation function in shear flow and obtained $`\tau _\alpha (T,\dot{\gamma })^1=\tau _\alpha (T,0)^1[1+A_\alpha \dot{\gamma }`$\] with $`A_\alpha \tau _\alpha (T,0)`$.

The sensitive shear-dependence of $`\tau _\alpha (T,\dot{\gamma })`$ predicted by (6.4) suggests potential importance of dielectric measurements in shear flow. As a first experiment, Matsuyama et al. measured the dielectric loss function $`ϵ^{\prime \prime }(\omega ,\dot{\gamma })`$ in steady shear $`\dot{\gamma }`$ in oligostyrene and polyisoprene melts. In the former melt at $`T=42^{}`$C, $`ϵ^{\prime \prime }(\omega ,\dot{\gamma })`$ decreased nonlinearly as a function of $`\dot{\gamma }`$ at low frequencies $`(\omega 10^5`$ s<sup>-1</sup>) in the shear-thinning regime $`(\dot{\gamma }15`$ s<sup>-1</sup>). Their finding indicates that $`\tau _\alpha `$ decreases as a function of $`\dot{\gamma }`$ in the non-Newtonian regime, consistently with (6.4), More systematic dielectric measurements in supercooled systems under shear flow are very informative.

## VII Summary

We have performed very long MD simulations of a supercooled polymer melt composed of $`M=100`$ short chains with bead number $`N=10`$ in quiescent and sheared conditions. We here summarize our main simulation results together with remarks.

(i) The stress relaxation function $`G(t)`$ is shown to follow a stretched exponential decay (1.1) on the scale of the $`\alpha `$ relaxation time $`\tau _\alpha `$ and then the Rouse relaxation (3.6) on the scale of $`\tau _\mathrm{R}`$.

(ii) The nonlinear shear regime sets in at extremely small shear rate of order $`\tau _\mathrm{R}^1`$ in supercooled states, where marked shear-thinning and shape changes of chains are found. Scattering and birefringence experiments from weakly sheared melts near $`T_\mathrm{g}`$ seem to be very promising.

(iii) In the nonlinear shear regime, each chain undergoes random tumbling in our melt as in the case of isolated polymer chains in shear flow. It is of great interest how this effect is universal in solutions and melts and how it influences macroscopic rheological properties. For example, we are interested in whether or not such tumbling occurs in sheared entangled polymers.

(iv)Transient stress divided by $`\dot{\gamma }`$ after application of shear flow obeys the linear growth $`_0^t𝑑t^{}G(t^{})G_0t`$ for strain less than 0.1 and then saturates into a non-Newtonian steady-state viscosity. This initial growth is much steeper than that predicted by the Rouse model. As a result, the stress-optical relation does not hold transiently under deformation even in the linear (zero-shear) limit. Its violation is more enhanced for larger shear rates. These are consistent with the experiments.

(v) The time-correlation functions in shear flow are calculated for the end-to-end vector and the modified particle displacement in (6.1). The former represents the relaxation of chain conformations, while the latter the monomeric relaxation on the spatial scale of the particle distance. We can then determine the shear-dependent relaxation times, $`\tau _\mathrm{R}(T,\dot{\gamma })`$ and $`\tau _\alpha (T,\dot{\gamma })`$. They decrease nonlinearly and behave differently as functions of $`\dot{\gamma }`$ in the nonlinear shear regime as in (6.3) and (6.4). It is then of great interest how these times behave in strong shear for much larger $`N`$. We conjecture that if $`N`$ is sufficiently large, shear should first influence the overall chain conformations, while it does not much affect the monomeric relaxations. We propose dielectric measurements in shear flow, which should give information of $`\tau _\alpha (T,\dot{\gamma })`$. In addition, as demonstrated in Fig.13, the effect of shear on the van Hove self-correlation function is analogous to raising the temperature above $`T_\mathrm{g}`$ as in supercooled binary mixtures.

## Acknowledgments

We would like to thank T. Inoue for valuable discussions on the stress-optical relation. This work is supported by Grants in Aid for Scientific Research from the Ministry of Education, Science, Sports and Culture of Japan. Computations have been performed at the Human Genome Center, Institute of Medical Science, University of Tokyo.


text:
# Observation of an incoherent thermally activated proton hopping process in calix-[4]-arene by means of anelastic spectroscopy.

## I INTRODUCTION

Calix\[n\]arenes are cup-like \[ln\] metacyclophanes which derive from the condensation of phenols and formaldehyde in different conditions. The bracketed number indicate the number of phenol units (the atoms inside the rectangle in Fig.1) and hence defines the size of the macrocycle. The tetramer is named calixarene (see Fig.1) The number of the phenolic units can vary from 3 to 14 and different substituents can be attached to the aromatic ring in the para position (position 2 in Fig.1) and at the phenolic oxygen. During the last decades calix\[n\]arene have been extensively utilised as molecular platforms for the synthesis of potent receptors for neutral and charged species and for their ability to form endo-cavity inclusion complexes with guests of complementary size. On the other hand the discovery of quantum tunnelling processes involving the hydrogen atoms of the OH groups was recently reported.  Both the design of molecules with an internal cavity able to recognize metal cations and/or neutral molecules  and the study of the transfer of protons involved in hydrogen bonding, which is fundamental to many chemical and biological processes , are important subjects; however, many questions remain unsolved and, especially concerning the second item, very few studies have been reported for the condensed state.

Quite recently NMR relaxometry measurements indicated that hydrogen atoms in calixarene and p-tert-butyl-calixarene might move both by means of a thermally activated hopping process and by a coordinated quantum tunneling process by which the four H atoms of the O-H groups move in a coordinated fashion giving rise to an interconversion between the two tautomers of the molecule.  In p-tert-butyl-calixarene the coordinated proton transfer (relaxation) rate , $`\tau _c`$, has been directly determined and it was found to be independent of temperature between 15 and 19 K, clearly demonstrating that the dynamical behavior of the system is non-Arrhenius in character.  In the calixarene molecule, in addition, a thermally activated proton hopping process was found, with a classical Arrhenius behavior described by a relaxation time $`\tau _c`$=1.2$``$10$`{}_{}{}^{12}`$e$`^{\frac{900K}{kT}}`$ s, which has been ascribed to rotational flips of the single O-H group (jump of the H atom between state 1 and 1’ in Fig.1).

In the following we will report anelastic spectroscopy measurements on calixarene containing four phenolic units (Fig.1). All the carbon and oxygen bonds are saturated by hydrogen atoms. We will show that a peak in the anelastic spectrum can be possibly ascribed to the same hopping process detected by NMR measurements and we will discuss in detail a quantitative analysis of this peak.

NMR relaxometry is a well-established experimental technique for the study of molecular dynamics. From the present study we can deduce that anelastic spectroscopy has great potentialities to study the movements of atoms in molecules at a mesoscopic (several-molecule) level.

## II EXPERIMENTAL AND RESULTS

Calixarene is produced as a powder sample. In order to measure its elastic properties, the powder was mixed with a pure KBr powder and was pressed at high pressure in order to obtain a bar of 7x38x2.0 mm3. A pure KBr bar was prepared in the same way and its anelastic spectrum was also measured for comparison.

The complex Young’s modulus $`E\left(\omega \right)=E^{}+iE^{\prime \prime }`$, whose reciprocal is the elastic compliance $`s=E^1`$, was measured as a function of temperature by electrostatically exciting the lowest odd flexural mode. The vibration amplitude was detected by a frequency modulation technique. The vibration frequency, $`\omega /2\pi `$, is proportional to $`\sqrt{E^{}}`$, while the elastic energy loss coefficient (or reciprocal of the mechanical $`Q`$) is given by $`Q^1(\omega ,T)=`$ $`E^{\prime \prime }/E^{}=`$ $`s^{\prime \prime }/s^{}`$, and was measured by the decay of the free oscillations or the width of the resonance peak. The imaginary part of the dynamic susceptibility $`s^{\prime \prime }`$ is related to the spectral density $`J_\epsilon (\omega ,T)=`$ $`𝑑te^{i\omega t}\epsilon \left(t\right)\epsilon \left(0\right)`$ of the macroscopic strain $`\epsilon `$ through the fluctuation-dissipation theorem, $`s^{\prime \prime }`$ $`\left(\omega /2k_\text{B}T\right)J_\epsilon `$.

In order to study the presence of thermally activated peaks in the anelastic spectrum of calixarene the measurements are carried out at different vibration frequencies. Higher order flexural modes could not be exited on the sample presently used; however, as the flexural mode frequency, $`\omega /2\pi `$, is directly proportional to the height of the sample, $`h`$, the specimen thickness was reduced with sand paper changing $`h`$ from 2.2 to 1.1 mm and consequently $`\omega /2\pi `$ changed by almost a factor 2.

The anelastic spectra of the pure KBr and the mixed KBr-calixarene bar are reported in Fig.2. In the temperature range from 1.1 to 400 K the KBr sample presents just a featureless background which increases slightly between 1 and 300 K and more significantly above room temperature. Instead, the calixarene sample exhibits three well developed peaks around 60 (peak P1), 150 (P2) and 270 K (P3). As they are not present in the KBr specimen, they can be attributed to intrinsic anelastic processes inside the calixarene molecule.

An analysis of the three peaks is presented in the next section, where it will be shown that the peak at the lowest temperature can be identified with the same thermally activated hopping process of H atoms detected by NMR measurements, giving rise to rotational flip of the OH bond.

## III DISCUSSION

In order to obtain useful information from the three peaks present in the spectrum of calixarene a quantitative analysis of them is needed. The most common model describing relaxation processes (anelastic, dielectric, etc.) is due to Debye. The contribution of a relaxation process to the dynamic Young’s modulus is:

$`\left[\begin{array}{c}Q^1\hfill \\ \delta E^{}/E^{}\hfill \end{array}\right]`$ $`=`$ $`Ev_0c(\lambda _1\lambda _2)^2`$ (6)
$`{\displaystyle \frac{f_1f_2}{T}}{\displaystyle \frac{1}{1+\left(\omega \tau \right)^{2\alpha }}}\left[\begin{array}{c}\left(\omega \tau \right)^\alpha \hfill \\ 1\hfill \end{array}\right]`$
where $`c`$ is the molar concentration of relaxing units; $`f_1`$ and $`f_2`$ are the equilibrium fractions in the two possible configurations $`1`$ and $`2`$ between which the relaxation can occur with a characteristic time $`\tau \left(T\right)`$ which classically follows an Arrhenius law, $`\tau =\tau _0\mathrm{exp}(E_a/kT)`$; $`\lambda _i`$ is the elastic dipole of the state $`i`$, i.e. the average macroscopic strain which the sample would have if each cell (molecule) contained one relaxing unit in the state $`i`$; $`v_0`$ is the cell (molecule) volume and $`\alpha 1`$. The contribution of a relaxation process to the absorption is a peak whose maximum is centered at the temperature at which $`\omega \tau =1`$. Since $`\tau `$ decreases as $`T`$ increases, the peak shifts to higher temperature if measured at higher $`\omega `$. Around the same temperature, the real part of the Young’s modulus presents a negative step whose amplitude, from theoretical models, is twice the peak height.

The Debye formula describes a relaxation process between two possible states which occurs with a single activation energy, $`E_a`$, and a single characteristic time, $`\tau `$. Often real processes are broader than Debye peaks; in those cases a correction factor, $`\alpha `$, is introduced in (1) (Fuoss-Kirkwood model ). Physically this means that the relaxation process has a distribution of activation energies and relaxation times which become broader as $`\alpha `$ decreases.

A first quantitative analysis of the $`Q^1`$ spectrum of calixarene was obtained after the subtraction of a smooth background, adding a Debye contribution for peak P1, a Debye contribution for peak P2 and a Fouss-Kirkwood peak for peak P3. The experimental data and the theoretical curves at the two measured frequencies are reported in Fig.3. The parameters of the three peaks are reported in Table I.

P1 and P2 can be described by the Debye process characterized by a single relaxation time, with the assumption that the $`\tau \left(T\right)`$ follows an Arrhenius law, $`\tau =\tau _0\mathrm{exp}(E_a/kT)`$. For peak P3 the fit is less satisfactory and this may be due to an imperfect substraction of the background, which in that range of $`T`$ starts to increase very fastly. However P3 does not seem to be larger than a pure Debye relaxation with the given parameters. For a Debye process, the height of the peak decreases as the frequency increases almost as 1/T, unless the atom or the defect relaxes between the states of an asymmetric potential well. In the present case for both higher temperature peaks the 1/T law is not valid and a multiplicative factor was introduced in (1) to account for the different height of the peaks at the two different frequencies. In the following we will not analyze the dependence of the peak heights on frequency, as the mechanical process used to obtain a second vibrational frequency could have affected the distribution of calixarene powders.

All the peaks are certainly due to some atomic motion inside the calixarene molecule as they are not present in the pure KBr specimen. From Table I one can notice that the $`\tau _0`$ values for the three peaks are the same within one order of magnitude. Peak P1 can be ascribed to a relaxation process involving the reorientation of a single bond or equivalently to the motion of a single atom. The similarity among the three characteristic times may suggest that also the higher temperature relaxation involves the motion of single atoms or a small number of them. Peak P3 is not present in the anelastic spectrum of p-tert-butyl-calixarene, whilst a peak similar to P2 is present also in the last compound. Considering the similarities of the structure of the calixarene and p-tert-butyl-calixarene molecules, it seems that peak P2 may be due to the reorientation of the H atoms at positions 2 and 2’ in Fig.1. However more experimental studies are needed to ascribe peak P2 and P3 to some physical processes. We want to point out that the activation energy and the relaxation time obtained from the fit procedure are not affected by the possible uncertainty on the peak height.

In the following we will deal with peak P1. It can be quite satisfactorily described by a Fuoss-Kirkwood model with $`\alpha =0.48`$; it means that the characteristic time and activation energy have a continuous distribution function around the values $`\tau _0=0.810^{12}s`$ and $`E_a=1120K`$. These values are quite similar to the parameters of the thermally activated proton hopping process, with a classical Arrhenius behavior measured by NMR relaxometry, which gave $`\tau _0=1.210^{12}s`$ and $`E_a=900K`$. This close similarity indicates that the NMR and anelastic relaxation processes may be originated by the same physical process, that is the rotational flips of the O-H groups of the calixarene molecules.

In order to gain insight and to have more information about the distribution functions of $`\tau `$ and $`E_a`$, we tried a different approach based on a direct evaluation of the distribution of activation energies; at this step we fixed $`\tau _0`$ to the single value obtained for the Fuoss-Kirkwood peak, both for simplicity of calculation and considering that we are dealing with a single relaxation process which corresponds to a well defined reorientation. This approach is very useful in describing a similar peak found in p-tert-butyl-calixarene, which cannot be modelled by a Debye or by a Fouss-Kirkwood relaxation or even by the sum of a big number of them. The expression for the elastic energy loss coefficient used for the fit procedure was:

$`Q^1`$ $`=`$ $`a{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{f_1f_2}{T}}`$ (8)
$`{\displaystyle \frac{\omega \tau \left(E_a\right)}{1+\left[\omega \tau \left(E_a\right)\right]^2}}g\left(E_a\right)dE_a`$
where $`a`$ is a constant, $`g\left(E_a\right)`$ is a normalized distribution function of the activation energy, $`E_a`$. Both a Lorentzian and a Gaussian distribution function have been tested. The Gaussian distribution is able to better reproduce the experimental data than a Lorentzian curve. The best fit curves using a Gaussian distribution function are reported in Fig.4. The best fit parameters are: $`\tau _0=0.810^{12}s`$, $`E_a=1070K`$, $`\sigma \left(E_a\right)=130K`$, where $`\sigma \left(E_a\right)`$ is the variance of the Gaussian distribution.

The existence of a distribution of activation energies more than a single $`E_a`$ can be explained considering that each H atom of the OH groups in a calixarene molecule moves in an environment (potential well) which can vary as the others three H atoms of the OH groups in the same molecule change their position. Moreover the potential energy of the H atom can be at least partially influenced by the surrounding calixarene molecules.

A distribution of activation energies for a relaxation process measured by NMR relaxometry, has been discussed also in the case of 1,3-di-t-butylbenzene, which was a very interesting case for the study of H dynamics in organic molecules in the condensed state.  The experimental data were described by a non-Bloembergen-Purcell-Pound (BPP) spectral density. In the case of anelastic relaxation, the BPP spectral density has its exact counterpart in the Debye model with an Arrhenius-like dependence of $`\tau `$ on $`T`$ and $`E_a`$. Beckmann et al. pointed out that the non-BPP behavior can also be interpreted in terms of a distribution of exponential correlation function (leading to the Arrhenius law), with the distribution at least partially characterized by a correlation time $`\tau `$ which could either be a cutoff correlation time or a mean correlation time, depending on the model. Correspondingly there is a distribution of activation energies. However the relation between the two distribution functions cannot be determined until the dependance of $`\tau `$ on $`E_a`$ is modeled. In the present case $`\tau _0`$, which has been fixed to the best fit value, can be interpreted as a mean correlation (relaxation) time. Moreover the Arrhenius dependence of $`\tau `$ on $`T`$ and $`E_a`$ has been explicitly used. However the Gaussian distribution function of the activation energies has not been derived from a model, but has been found to describe quite well the experimental data.

A last consideration regards the lack of any quantum tunneling process in the anelastic spectrum of calixarene, as that detected by NMR relaxometry. One possible explanation for this, is that we perform measurements of the dynamical Young’s modulus as a function of temperature, measured with a practically fixed frequency. If the correlation time $`\tau `$ in the tunnelling process becomes temperature independent, like in the case of p-tert-butyl-calixarene, the condition $`\omega \tau =1`$ in order to have a peak in the elastic energy loss coefficient could be never fulfilled.

## IV CONCLUSION

We measured the anelastic spectrum of a bar of KBr containing a few percent of calixarene powder. We observed three thermally activated peaks. The lowest temperature peak can be described both by a Fuoss-Kirkwood model and by a model considering a Gaussian distribution of activation energies, leaving the correlation time $`\tau _0`$ fixed. Due to the similarity between the physical parameters of this peak and of a relaxation found in the NMR spectrum of the same compound, the anelastic peak can be ascribed to hopping of H atoms of the OH groups, which corresponds to a flip-flop of the OH bond. From the results of the present study, anelastic spectroscopy appears as a good tool for studying atomic motion inside molecules at a mesoscopic (few molecules) level.


text:
# Lipid membranes with an edge

## I Introduction

A closed lipid membrane is described remarkably well by a geometrical hamiltonian. This hamiltonian is constructed as a sum of the scalars, truncated at an appropriate order, which characterize those features of the membrane geometry which are relevant. A term quadratic in the extrinsic curvature provides a measure of the energy penalty associated with bending Nelson et al. (1989); Peliti (1996); Safran (1994); Lipowsky and Sackmann (1995); Seifert (1997); any intrinsic tendency to bend one way and not the other is captured by a term linear in the extrinsic curvature Helfrich (1973).

The shape equation determining the equilibria of this membrane is a fourth-order non-linear elliptic PDE of the form $`^2K+`$$`K^3`$$`=0`$, where $`^2`$ is the Laplacian on the membrane, $`K`$ is the sum of the principal curvatures, and by $`\mathrm{`}K^3`$’ we mean a cubic polynomial in these curvatures Zhong-Can and Helfrich (1987). Here, we would like to examine the boundary conditions which must supplement this equation when the membrane possesses a free edge. The energy cost associated with this edge is, to a good approximation, proportional to the exposed length. During the formation process, material will either be added to the edge or the edge will heal itself so as to form closed structures. However, unlike a simple soap film with a free edge, a lipid membrane does possess stable (cup shaped) equilibria Boal and Rao (1992). An examination of such structures is a small but important step towards a fuller understanding of self assembly.

Our primary focus will be on the boundary geometry. We have a surface with a boundary and a certain energy penalty associated with it, a well defined problem in classical field theory. The boundary conditions are identified by demanding that the energy be stationary for arbitrary deformations of the edge geometry. In distinction to earlier analytic treatments of this problem, we will relax the assumption that the membrane geometry be axially symmetric Boal and Rao (1992). This is important not only for conceptual reasons. Generally, there will be no priveleged parametrization such as that tailored to axial symmetry; in an axially symmetric geometry the edge itself is simply a circle. We find that there are three boundary conditions. As we will demonstrate one of these conditions, involving three derivatives of the embedding function, is satisfied identically in the axially symmetric limit. Therefore this limit cannot be considered as a reliable guide to the general case.

While the variational approach does capture the geometrical nature of the boundary conditions, the physical interpretation of these conditions still needs to be clarified. Ideally, one would like to interpret them in terms of the balance of forces operating at the edge. To do this in a way which does justice to the geometry, we identify the conserved Noether currents associated with the intrinsic translational invariance of the configuration Capovilla and Guven (2002a). The three, apparently unrelated, boundary conditions are now cast in terms of the three components of a single vector identity on the edge.

We finish with a discussion of the effect of a Gaussian rigidity term on a lipid membrane with edges. Whereas such a topological term does not alter the bulk shape equation, we show that it does modify the boundary conditions which apply to it in a way which will have consequences in the bulk. This extension is relevant in topology changing processes.

The outline of the paper is as follows. In Sect. II we consider the simple example of a surface tension dominated membrane. This allows us to establish our notation and to derive the boundary conditions in a simple context. In Sect. III, we derive the boundary conditions at the edge for a lipid membrane. We then specialize to axially symmetric configurations to compare our results in this limit with previous work on the subject. In Sect. IV, we consider the balance of the forces operating at the edge and we show how they are related to the boundary conditions. The effect of a Gaussian rigidity term to the membrane energy is the subject of Sect. V, where we obtain the appropriate modifications in the boundary conditions. We end with some remarks in Sect. VI.

## II Surface tension pitted against edge tension

It is useful to examine first the simpler situation in which the membrane physics is dominated by surface tension, such as a soap film with a free edge. Let the membrane surface $`\mathrm{\Sigma }`$ have an area $`A`$, with boundary $`𝒞`$ of length $`L`$, and the tension in the membrane bulk be a constant $`\mu `$, and that on the edge $`\sigma `$. The energy is then given as a sum of two terms

$$F=\mu A+\sigma L.$$
(1)
Surface tension tends to decrease the membrane area; line tension to decrease the length of the free boundary. Without some further refinement this model does not admit stable equilibria. Suppose a hole is punctured in the film, then depending on its radius, either the hole will close healing the film, or grow and destroy it. An unstable equilibrium clearly exist when the radius is tuned to coincide with a critical value $`r_c`$. On dimensional ground, one world expect $`r_c\mu /\sigma `$. (We will ignore stability here as our interest in this model is only as a point of reference for a lipid membrane.)

The membrane surface $`\mathrm{\Sigma }`$ is described by the embedding $`𝐗`$ in three dimensional space $`R^3`$ as $`𝐱=𝐗(\xi ^a)`$, where $`𝐱`$ are coordinates for $`R^3`$, and $`\xi ^a`$ coordinates for the surface ($`a,b,\mathrm{}=1,2`$). Its edge $`𝒞`$ is embedded in turn as a curve on $`\mathrm{\Sigma }`$ as $`\xi ^a=Y^a(s)`$ which we parametrize by its arclength $`s`$. We can now cast $`F`$ as

$$F=\mu _\mathrm{\Sigma }d^2\xi \sqrt{\gamma }+\sigma _𝒞𝑑s.$$
(2)
Here, the metric induced on $`\mathrm{\Sigma }`$ is given by $`\gamma _{ab}=𝐞_a𝐞_b`$, where $`𝐞_a:=_a𝐗`$ are tangent to the surface, $`\gamma =\mathrm{det}\gamma _{ab}`$, and $`dA=\sqrt{\gamma }d^2\xi `$. Note that we can also consider the direct embedding of the edge $`𝒞`$ in $`R^3`$, via $`𝐱=𝐘(s)`$, where $`𝐘=𝐗(Y^a(s))`$. The tangent to $`𝒞`$ in $`R^3`$ is equivalently expressed in either of two ways: $`𝐭=𝐞_at^a`$, where $`t^a=\dot{Y}^a`$; or $`𝐭=\dot{𝐘}`$, where a dot denotes a derivative with respect to arclength $`s`$.

The energy is a functional of the embedding $`𝐗`$ of $`\mathrm{\Sigma }`$ in $`R^3`$. There is no need to vary the edge embedding $`Y^a`$ independently: the $`Y^a`$ are fixed by the constraint that the two embeddings agree, $`𝐘=𝐗`$, on $`𝒞`$.

Equilibrium configurations are those at which the energy Eq.(1) is stationary. To derive the equations describing the equilibrium configurations of a membrane described by this model, we first consider a variation of the embedding $`𝐗`$ of the membrane $`𝐗𝐗+\delta 𝐗`$. We let $`𝐧`$ denote the unit normal to the surface $`\mathrm{\Sigma }`$. We decompose the displacement with respect to the spatial basis adapted to $`\mathrm{\Sigma }`$, $`\{𝐞_a,𝐧\}`$, as,

$$\delta 𝐗=\mathrm{\Phi }^a𝐞_a+\mathrm{\Phi }𝐧.$$
(3)
We now have that the induced metric varies according to (see e.g. Capovilla and Guven (2002a))

$$\delta _X\gamma _{ab}=2K_{ab}\mathrm{\Phi }+_a\mathrm{\Phi }_b+_b\mathrm{\Phi }_a,$$
(4)
where $`K_{ab}`$ denotes the extrinsic curvature tensor,

$$K_{ab}=𝐞_b_a𝐧,$$
(5)
and $`_a`$ is the covariant derivative on $`\mathrm{\Sigma }`$ compatible with $`\gamma _{ab}`$. The derivative terms in the variation of $`\gamma _{ab}`$ are its Lie derivative along the tangential vector field, $`\mathrm{\Phi }^a`$. The variation of $`A`$ is

$`\delta _XA`$ $`=`$ $`{\displaystyle _\mathrm{\Sigma }}𝑑A\left[K\mathrm{\Phi }+_a\mathrm{\Phi }^a\right]`$ (6)
$`=`$ $`{\displaystyle _\mathrm{\Sigma }}𝑑AK\mathrm{\Phi }+{\displaystyle _𝒞}𝑑sl^a\mathrm{\Phi }_a.`$
The mean extrinsic curvature, $`K=K_{ab}\gamma ^{ab}`$. The second line follows from the preceding one using Stoke’s theorem. Here $`l^a`$ is the outward pointing normal to $`𝒞`$ on $`\mathrm{\Sigma }`$. Only the normal projection $`\mathrm{\Phi }`$ of the variation plays a role in determining the bulk equilibrium of the membrane. This is true generally, regardless of the model. In this particular model, however, there is no boundary term associated with the bulk normal displacement $`\mathrm{\Phi }`$. As we will see, this is not generally true — a happy accident when the energy is truncated at the area term. On the other hand, the tangential bulk variation $`\mathrm{\Phi }^a`$ always gives only a boundary term. This is a consequence of the fact that a tangential deformation corresponds in the bulk to an infinitesimal reparametrization of the surface. There is however a physical displacement of the boundary. In fact, the boundary contribution to Eq.(6) is easily identified as the change in the surface area of $`\mathrm{\Sigma }`$ under a normal deformation of its boundary, $`\delta Y^a=(l^b\mathrm{\Phi }_b)l^a`$, which at each point is directed along the tangent plane of $`\mathrm{\Sigma }`$ at that point. The projection of $`\mathrm{\Phi }^b`$ onto the edge $`𝒞`$ itself, $`t^a\mathrm{\Phi }_a`$, does not contribute.

Let us turn now to the variation of the edge embedding $`𝐘`$ induced by the bulk variation $`\delta 𝐗`$. It can be decomposed with respect to a basis adapted to both embeddings, $`𝐗`$ and $`𝐘`$, given by $`\{𝐭,𝐥,𝐧\}`$, where $`𝐥=𝐞_al^a`$. Thus at the edge we set

$$\delta 𝐘=\varphi 𝐭+\psi 𝐥+\mathrm{\Phi }𝐧,$$
(7)
where the edge and bulk components are identified by continuity, $`\psi =l^a\mathrm{\Phi }_a`$ and $`\varphi =t^a\mathrm{\Phi }_a`$. Modulo a divergence associated with a reparametrization of the boundary which involves the tangential component $`\varphi `$ that we can safely discard, we have for the variation of the infinitesimal arclength

$$\delta _Yds=ds\left(\kappa \psi +K_{}\mathrm{\Phi }\right),$$
(8)
where we have used the fact that

$$\dot{𝐭}=\kappa 𝐥K_{}𝐧.$$
(9)
Here $`\kappa `$ is the geodesic curvature of $`𝒞`$ associated with its embedding in $`\mathrm{\Sigma }`$, and we have defined $`K_{}=K_{ab}t^at^b`$. The unconventional minus sign in the first term of Eq. (9) comes about because $`𝐥`$ is the outward normal to $`𝒞`$ on $`\mathrm{\Sigma }`$, i.e. $`\dot{t}^a=\kappa l^a`$.

The corresponding deformation in $`L`$ is then given by

$$\delta _YL=_𝒞𝑑s\left(\kappa \psi +K_{}\mathrm{\Phi }\right).$$
(10)
Summing the two contributions (6) and (10) to the variation of the energy $`F`$, as given by (1), we find

$$\delta _XF=\mu _\mathrm{\Sigma }𝑑AK\mathrm{\Phi }+_𝒞𝑑s\left[(\mu +\sigma \kappa )\psi +\sigma K_{}\mathrm{\Phi }\right].$$
(11)
The bulk equilibrium is a minimal surface unaffected by the boundary, satisfying $`K=0`$. On the boundary, the projections along the normals to the edge, $`\psi `$ and $`\mathrm{\Phi }`$, represent independent deformations, so that stationarity of $`F`$ requires the vanishing of the corresponding coefficents. We thus read off the two boundary conditions:

$`\sigma \kappa +\mu `$ $`=`$ $`0,`$ (12)
$`\sigma K_{}`$ $`=`$ $`0.`$ (13)
The first tells us that the geodesic curvature of the edge as embedded in the membrane is constant. The second simply enforces the vanishing of $`K_{}`$ at the edge. Note that the completeness of the basis $`\{𝐭,𝐥\}`$ of tangent vectors on $`\mathrm{\Sigma }`$ at $`𝒞`$, $`\gamma ^{ab}=t^at^b+l^al^b`$, permits us to express the mean curvature at the edge as $`K=K_{}+K_{}`$, where $`K_{}=K_{ab}l^al^b`$. Thus modulo the bulk equilibrium $`K=0`$, the boundary condition (13) can be alternatively expressed as $`K_{}=0`$. The only potentially non vanishing component of $`K_{ab}`$ on the edge is the off-diagonal component, $`K_{}=t^al^bK_{ab}`$.

For this particular model our approach has been heavyhanded; the boundary conditions we have written down are an elaborate way to express the simple vector identity

$$\sigma \dot{𝐭}=\mu 𝐥,$$
(14)
which equates the change in the tension over the interval $`\mathrm{\Delta }s`$ along the edge, $`\sigma \mathrm{\Delta }𝐭`$, to the force due to surface tension acting on the edge, $`\mu 𝐥\mathrm{\Delta }s`$. The apparent mismatch in counting (three versus two) is accounted for by noting that the projection of Eq.(14) along $`𝐭`$ is an identity. For higher order models, as we will see, this projection will not be vacuous.

Note that had we $`N`$ sheets conjoined on a single edge, Eq.(14) gets modified in an obvious way:

$$\sigma \dot{𝐭}=\mu \underset{i=1}{\overset{N}{}}𝐥_i,$$
(15)
where $`𝐥_i`$ is the vector normal to the edge which is tangent to the $`i^{\mathrm{th}}`$ sheet. Eq.(15) provides a generalization of the Neumann rule for soap bubble clusters at a Plateau border Weaire and Hutzler (1999) to accomodate line tension on the edge. A simple application is considered in Capovilla and Guven (2002b).

## III Lipid membrane with an edge

A lipid membrane is modeled by a phenomenological energy quadratic in the extrinsic curvature of the surface. Let us write this as

$$F_b=_\mathrm{\Sigma }𝑑A(\gamma ^{ab},K_{ab}),$$
(16)
i.e. $``$ depends at most on the extrinsic curvature, and not, for example, on its derivative $`_aK_{bc}`$. In particular, we will focus on the model described by the Helfrich energy density

$$=\alpha (KK_0)^2+\mu .$$
(17)
The spontaneous curvature $`K_0`$ is a constant, as is the bending rigidity $`\alpha `$. The constant $`\mu `$ is interpreted here as the Lagrange multiplier implementing the constraint on the membrane area. We will discuss the addition of a Gaussian rigidity term in the next section.

The energy of the bulk and the edge is

$$F=F_b+\sigma L.$$
(18)
The shape equation describing the equilibrium in the bulk, which is derived from the extremization of the energy (17)

$$\alpha [2^2K+2(KK_0)+(K_0{}_{}{}^{2}K^2)K]+\mu K=0,$$
(19)
is well known Zhong-Can and Helfrich (1987). The structure of this equation has been discussed in detail elsewhere Capovilla and Guven (2002a), where a novel derivation is also provided. The scalar curvature $``$ appearing in Eq.(19) is related to the extrinsic curvature through the Gauss-Codazzi equation,

$$=K^2K_{ab}K^{ab}.$$
(20)
Under a tangential deformation of the surface, $`\delta _{}𝐗=\mathrm{\Phi }^a𝐞_a`$, the energy density transforms as a divergence which is transferred to the boundary,

$$\delta _{}F_b=_𝒞𝑑sl_a\mathrm{\Phi }^a.$$
(21)
This is because the local scalar energy density $``$ transforms as

$$\delta _{}=\mathrm{\Phi }^a_a.$$
(22)
The details of $``$ do not enter. Note that Eq.(21) agrees with the corresponding expression for the area with $`=1`$. As before, this boundary term induces a source into the boundary Euler-Lagrange equation. For an edge with a line tension $`\sigma `$, we get the first boundary condition, due to a deformation along the normal $`𝐥`$, $`\psi `$,

$$\sigma \kappa +=0,$$
(23)
where we have used Eqs. (10) and (21). This should be compared with Eq.(12) to which it reduces if $`=\mu `$, a constant. This boundary condition relates the geometry of the edge to the extrinsic curvature of the membrane evaluated at the boundary.

We now examine a normal deformation of the surface $`\mathrm{\Sigma }`$, $`\delta _{}𝐗=\mathrm{\Phi }𝐧`$. The shape equation (19) determining the local membrane equilibrium is obtained by demanding that the energy be stationary with respect to normal deformations of $`\mathrm{\Sigma }`$ which may or may not vanish on the boundary. As such this equation cannot be affected by the addition of a boundary. To determine the boundary conditions we need to extend the support of the variation to include the boundary. We have that the normal variation of the bulk energy can be written as

$$\delta _{}F_b=𝑑A\left[K\mathrm{\Phi }+2\mathrm{\Gamma }_{ab}K^{ab}\mathrm{\Phi }+^{ab}\delta _{}K_{ab}\right],$$
(24)
where $`\mathrm{\Gamma }_{ab}=/\gamma ^{ab}`$ and $`^{ab}=/K_{ab}`$. The boundary term we wish to identify in $`\delta _{}F_b`$ originates in the $`\delta _{}K_{ab}`$ term in this expression. We recall that the extrinsic curvature transforms as follows under a normal deformation of $`\mathrm{\Sigma }`$ (see e.g. Capovilla and Guven (2002a)):

$$\delta _{}K_{ab}=_a_b\mathrm{\Phi }+K_{ac}K^c{}_{b}{}^{}\mathrm{\Phi }.$$
(25)
We thus have that

$$\delta _{}F_b=𝑑A\left[\mathrm{\Phi }+_a\left(\mathrm{\Phi }_b^{ab}^{ab}_b\mathrm{\Phi }\right)\right],$$
(26)
where we have defined the Euler-Lagrange derivative

$$=(_a_b+K_{ac}K^c{}_{b}{}^{})^{ab}+K+2\mathrm{\Gamma }_{ab}K^{ab}.$$
(27)
Thus, modulo the bulk shape equation, $`=0`$, we get that the boundary contribution is

$$\delta _{}F_b=_𝒞𝑑sl_a\left[\mathrm{\Phi }_b^{ab}^{ab}_b\mathrm{\Phi }\right].$$
(28)
The terms proportional to $`_a\mathrm{\Phi }`$ and $`\mathrm{\Phi }`$ are not independent: the projection of $`_a\mathrm{\Phi }`$ along the edge is completely determined once $`\mathrm{\Phi }`$ is specified on $`𝒞`$. To decompose $`\delta _{}F`$ into two independent parts we proceed as follows: we first decompose $`_a\mathrm{\Phi }`$ into its normal and tangential parts with respect to $`𝒞`$,

$$_a\mathrm{\Phi }=l_a_{}\mathrm{\Phi }+t_a\dot{\mathrm{\Phi }},$$
(29)
where we have defined $`_{}=l^a_a`$. We now perform an integration by parts on the $`\dot{\mathrm{\Phi }}`$ term to obtain for the second term on the right hand side of Eq.(28)

$`{\displaystyle _𝒞}`$ $`ds`$ $`l_a^{ab}_b\mathrm{\Phi }=`$
$`{\displaystyle _𝒞}`$ $`ds`$ $`\left[l_al_b^{ab}_{}\mathrm{\Phi }\mathrm{\Phi }{\displaystyle \frac{d}{ds}}\left(l_a^{ab}t_b\right)\right],`$ (30)
where we have discarded a total derivative term with respect to arclength. In this way we succeed in isolating the independent normal variations at the boundary, the coefficients of $`\mathrm{\Phi }`$ and $`_{}\mathrm{\Phi }`$.

From Eqs. (10), (28), (30), we obtain that the total boundary contribution of the normal variations is

$`\delta _{}`$ $`F_b`$ $`={\displaystyle _𝒞}ds\{l_al_b^{ab}_{}\mathrm{\Phi }`$
$`+`$ $`[`$ $`l_a_b^{ab}+{\displaystyle \frac{d}{ds}}\left(l_a^{ab}t_b\right)+\sigma K_{}]\mathrm{\Phi }\},`$ (31)
so that we can immediately read off the two boundary conditions which supplement Eq. (23),

$`l_a_b^{ab}+{\displaystyle \frac{d}{ds}}\left(l_a^{ab}t_b\right)+\sigma K_{}`$ $`=`$ $`0,`$ (32)
$`l_al_b^{ab}`$ $`=`$ $`0.`$ (33)
The first is of third order in derivatives of the embedding functions. This is consistent with the fact that the shape equation (19) is of fourth order. Using the decomposition of the covariant derivative (29), it can be written in the alternative form

$$l_al_b_{}^{ab}+2\frac{d}{ds}\left(l_a^{ab}t_b\right)+\kappa (l_al_bt_at_b)^{ab}+\sigma K_{}=0.$$
(34)
In the case of a membrane described by the Helfrich Hamiltonian (17) with an edge the third boundary condition (33) implies

$$K=K_0,$$
(35)
on the edge — the rigid membrane necessarily has a constant mean curvature at the edge equal to its spontaneous value. This is entirely independent of the tensions, $`\mu `$ or $`\sigma `$, or of the rigidity modulus. If $`K_0=0`$, the membrane is minimal at its edge.

The second boundary condition (32) is of Robin type. For any $``$ which is a function only of $`K`$, we have that $`^{ab}\gamma ^{ab}`$, so that the middle term in (32) vanishes,

$$l_a^{ab}t_b=0,$$
(36)
and the boundary condition reduces to

$$2\alpha _{}K+\sigma K_{}=0.$$
(37)
This equation determines the normal derivative of $`K`$ in terms of the component of the extrinsic curvature tangent to the edge. It does not involve the surface tension $`\mu `$. We emphasize that its existence seems to have gone unnoticed.

The first boundary condition, Eq.(23), together with Eq.(35), implies that on the edge

$$\sigma \kappa +\mu =0.$$
(38)
The geodesic curvature of a loaded boundary is completely fixed by the ratio of the tensions in exactly the same way as in the previous section for soap bubbles, see Eq. (12) . If $`\mu =0`$ the edge is necessarily a geodesic of the bulk geometry.

If the line tension on the boundary vanishes, $`\sigma =0`$, the consistency of Eq.(35) with (23) requires that $`\mu =0`$ also. Furthermore Eq.(37) implies $`_{}K=0`$ on the boundary. But the unique solution satisfying the two boundary conditions $`K=K_0`$ and $`_{}K=0`$ is $`K=K_0`$ everywhere. One way to see this is to construct the Gaussian normal coordinates adapted to the edge, $`(l,s)`$, where $`l`$ is the length of the geodesic which intersects the edge normally. With respect to this system of coordinates, the Laplacian assumes the form $`^2=_l^2+\kappa _l+_s^2`$ in the neighborhood of the edge. Thus, modulo the boundary conditions, $`^2K=_l^2K`$ on the edge. But Eq.(19) implies that $`^2K=0`$ there so that $`_l^2K`$ and all higher derivatives vanish. If $`K`$ is analytic in $`l`$, then $`K=K_0`$. If $`\mu 0`$, there is no such constraint. The geometry is very severely constrained by the boundary conditions.

Let us now examine an axially symmetric membrane with an axially symmetric edge. With respect to cylindrical polar cooordinates $`\{\rho ,z,\phi \}`$ on $`R^3`$, the membrane is described by $`\rho =R(\mathrm{})`$, $`z=𝒵(\mathrm{})`$ where $`𝒵^{}{}_{}{}^{2}+R^{}{}_{}{}^{2}=1`$. $`\mathrm{}`$ is the arc length along a a curve with fixed $`\phi `$, and the primes denote a derivative with respect to $`\mathrm{}`$. The intrinsic geometry of $`\mathrm{\Sigma }`$ is described by the line element

$$d\tau ^2=d\mathrm{}^2+R^2(\mathrm{})d\phi ^2,$$
(39)
We can write the extrinsic curvature in a form consistent with axial symmetry as

$$K_{ab}=\mathrm{}_a\mathrm{}_bK_{\mathrm{}}+(\gamma _{ab}\mathrm{}_a\mathrm{}_b)K_R,$$
(40)
where $`K_l`$ and $`K_R`$ are two spatial scalars which we identify as the principal curvatures of the embedding of $`\mathrm{\Sigma }`$ in $`R^3`$, $`\mathrm{}^a`$ is the outward pointing unit normal to the circle of fixed $`\mathrm{}`$, $`\mathrm{}^a=(1,0)`$. The mean curvature is $`K=K_{\mathrm{}}+K_R`$. To evaluate the principal curvatures, it is convenient to define $`\mathrm{\Theta }`$ as the angle which the tangent to a curve of fixed $`\phi `$ makes with the cylindrical radial direction:

$$\frac{d𝒵}{dR}=\mathrm{tan}\mathrm{\Theta }.$$
(41)
We then have $`𝒵^{}=\mathrm{sin}\mathrm{\Theta }`$, and $`R^{}=\mathrm{cos}\mathrm{\Theta }`$, so that the principal curvatures are

$$K_{\mathrm{}}=\mathrm{\Theta }^{},K_R=\frac{\mathrm{sin}\mathrm{\Theta }}{R}.$$
(42)
Axial symmetry implies that the fourth order shape equation can be integrated to provide a third order equation for $`R`$ as a function of $`\mathrm{}`$. It has been shown elsewhere (Zheng and Liu (1993); Boz̆ic̆ et al. (1997), see also Capovilla and Guven (2002a)) that this equation takes the form

$$2\alpha \mathrm{cos}\mathrm{\Theta }\left(\mathrm{\Theta }^{}+\frac{\mathrm{sin}\mathrm{\Theta }}{R}\right)^{}+\alpha (\mathrm{\Theta }^{}+\frac{\mathrm{sin}\mathrm{\Theta }}{R})(\mathrm{\Theta }^{}\frac{\mathrm{sin}\mathrm{\Theta }}{R})\mathrm{sin}\mathrm{\Theta }+2\alpha K_0\frac{\mathrm{sin}^2\mathrm{\Theta }}{R}(\mu +\alpha K_0^2)\mathrm{sin}\mathrm{\Theta }=0.$$
(43)
If the boundary $`𝒞`$ is also axially symmetric so that it coincides with a fixed value of $`\mathrm{}`$ then $`l^a=\mathrm{}^a`$, $`K_{}=K_R`$, $`K_{}=K_{\mathrm{}}`$ and $`K_{}=0`$. It is simple to show that $`\kappa =R^{}/R`$. We thus have for the boundary conditions, Eqs.(35) and (38),

$$\mathrm{\Theta }^{}+\mathrm{sin}\mathrm{\Theta }/R=K_0,\sigma R^{}=\mu R.$$
(44)
The remaining boundary condition, Eq.(37), of third order in derivatives appears to present a problem: a third order ODE does not admit third order boundary conditions. The inconsistency, however, is only apparent: on the boundary, the shape equation Eq.(43) itself reproduces, modulo Eqs.(44), the troublesome boundary condition Eq.(37). Our analysis is thus completely consistent with the axially symmetric analysis of Boal and Rao (1992) where the boundary conditions Eqs.(23) and (33) are derived. It is worth stressing, however, that potential pitfalls of using the axially symmetric problem as a guide to the more general problem. The boundary condition (32) is a non-trivial constraint on the geometry which is not already encoded in the shape equation for non-axially symmetric configurations.

## IV Balance of forces

In this section, we consider the balance of the forces operating at the edge. This provides the missing intuition on the physical origin of the boundary conditions we have derived in the previous section.

Consider a point on the edge. In equilibrium, the tension $`𝐠`$ must satisfy

$$\dot{𝐠}=𝐟^al_a.$$
(45)
Here $`𝐟^a`$ is the membrane stress tensor so that $`𝐟^al_a`$ is the ‘surface tension’ acting on the edge due to unbalanced stresses in the bulk at its boundary. In Capovilla and Guven (2002a), it was shown that the bulk stress tensor for the model defined by the Helfrich energy (17) can be expressed in the form

$`𝐟^a`$ $`=`$ $`[2\alpha K(K^{ab}{\displaystyle \frac{K}{2}}\gamma ^{ab})2\alpha K_0(K^{ab}K\gamma ^{ab})`$ (46)
$``$ $`(\mu +\alpha K_0^2)\gamma ^{ab}]𝐞_b2\alpha ^aK𝐧.`$
Thus its projection along the normal to the edge $`l^a`$ is

$`𝐟^al_a`$ $`=`$ $`\left\{2\alpha (KK_0)K_{}\alpha (KK_0)^2\mu \right\}𝐥`$ (47)
$`+`$ $`2\alpha (KK_0)K_{}𝐭2\alpha _{}K𝐧.`$
In addition, as we have seen in Sect. II,

$$𝐠=\sigma 𝐭.$$
(48)
Using Eq. (9) for $`\dot{𝐭}`$, we read off the three components of Eq.(45),

$`\sigma \kappa `$ $`=`$ $`2\alpha (KK_0)K_{}\alpha (KK_0)^2\mu ,`$ (49)
$`\sigma K_{}`$ $`=`$ $`2\alpha _{}K,`$ (50)
$`0`$ $`=`$ $`2\alpha (KK_0)K_{},`$ (51)
respectively along $`𝐥`$, $`𝐧`$, and $`𝐭`$. The condition (50) coincides with the boundary condition (37). If $`K_{}0`$, Eq.(51) implies that $`K=K_0`$. The remaining boundary condition (49) then coincides with a linear combination of the boundary conditions (35) and (38). In the axially symmetric geometry, however, $`K_{}`$ does vanish so that (51) does not imply $`K=K_0`$ as it stands. One needs then to appeal to the integrated shape equation (43), which together with (49) and (50) reproduce $`K=K_0`$.

We thus have identified a very simple (if heavily disguised) physical interpretation of the boundary conditions. In particular, in this approach, the boundary condition $`K=K_0`$ emerges as the vanishing of the stress induced by the bulk along the edge. Note that the variational approach did not rely on the identification of projections. Indeed, the boundary condition corresponding to the projection along $`𝐭`$ was originally identified by demanding stationary energy for independent boundary variations of $`_{}\mathrm{\Phi }`$.

We also note that the form of Eq.(45) implies the integrability condition,

$$_𝒞𝑑s𝐟^al_a=0$$
(52)
on the edge. The existence of these three extremely non-trivial conditions is far from obvious in our previous approach.

One can say more. Take the equation, $`_a𝐟^a=0`$, describing the conservation of the stress tensor. dot it into $`𝐗`$ and integrate over the membrane surface. We get

$$𝑑A_a(𝐗𝐟^a)=𝑑A𝐞_a𝐟^a.$$
(53)
Working on the right, we have

$`{\displaystyle 𝑑A_a(𝐗𝐟^a)}`$ $`=`$ $`{\displaystyle 𝑑s𝐗l_a𝐟^a}`$ (54)
$`=`$ $`{\displaystyle 𝑑s𝐗\dot{𝐠}}`$
$`=`$ $`{\displaystyle 𝑑s𝐭𝐠}`$
$`=`$ $`\sigma L,`$
where we have used Eq.(45) on the second line, as well as Eq.(48) on the last line. On the other hand, if we write $`𝐟^a=f^{ab}𝐞_b+f^a𝐧`$, then

$$dA𝐞_a𝐟^a=dAf^a{}_{a}{}^{}.$$
(55)
We note that the bending energy $`𝑑AK^2`$ is a conformal invariant, and so does not contribute to the trace $`f^a_a`$. We have $`f^a{}_{a}{}^{}=2(\alpha K_0(KK_0)\mu )`$, so that

$$2\mu A2\alpha K_0𝑑A(KK_0)+\sigma L=0.$$
(56)
This condition is useful for identifying the sign associated with the mulupliers. For example, if $`K_0=0`$, it is clear that $`\mu `$ is necessarily negative.

## V Gaussian rigidity

The geometrical scalars we can construct with dimension $`[\mathrm{Length}]^2`$ are $``$, $`K^2`$ and $`K_{ab}K^{ab}`$. The Gauss-Codazzi equation (20) tells us that the three scalar invariants $``$, $`K^2`$ and $`K_{ab}K^{ab}`$ are not independent. In addition, the Gauss-Bonnet functional

$$I=_\mathrm{\Sigma }𝑑A$$
(57)
is a topological invariant if the membrane is closed. More generally for an open membrane,

$$I=_\mathrm{\Sigma }𝑑A+2_𝒞𝑑s\kappa $$
(58)
is a topological invariant. A consequence is that if a Gaussian rigidity term is included in the energy a line rigidity $`𝑑s\kappa `$ is necessarily induced along its boundary.

To obtain the variation of the Gaussian term, we need to know how the scalar curvature varies. Its tangential deformation is straightforward; to determine its normal deformation, we exploit Eq.(20) and with it the technology developed Sects. II and III.

Consider now a Gaussian rigidity addition to the bulk energy, so that

$$F=F_b+\beta 𝑑A,$$
(59)
Whereas the bulk shape equation is unmodified, all three boundary conditions are changed:

$`\alpha (KK_0)^2+\beta +\mu +\sigma \kappa `$ $`=`$ $`0,`$ (60)
$`2\alpha _{}K2\beta \dot{K}_{}+\sigma K_{}`$ $`=`$ $`0,`$ (61)
$`\alpha (KK_0)+\beta K_{}`$ $`=`$ $`0.`$ (62)
We note that the Gauss-Codazzi equation (20) allows us to express $``$ in terms of the projections of $`K_{ab}`$ with respect to the edge, $`=2(K_{}K_{}K_{}{}_{}{}^{2})`$. Eq.(60) is quadratic in the extrinsic curvature. Eqs.(61) and (62) by constrast are linear relationship between $`K_{}`$ and $`K_{}`$.

Note that, unlike the case of the pure Helfrich model, the central term in Eq.(61) does not vanish in general. However, it does vanish in an axially symmetric geometry (with axially symmetric edge), $`K_{}=0`$. More generally, we have the integral statement

$$_𝒞𝑑s\left[2\alpha _{}K+\sigma K_{}\right]=0.$$
(63)
In an axially symmetric geometry one can check that, modulo the lower order boundary conditions (60) and (62), Eq.(43) reproduces Eq.(61) on the boundary.

Let us consider now the balance of the forces in this case. The Gaussian term makes no contribution to $`𝐟^a`$ Capovilla and Guven (2002a). Naively reinvoking Eq.(45) would appear to suggest that this term cannot modify the boundary conditions, in contradiction with what we have just derived. However, with a general function $``$ and in particular for $`\beta `$, the term $`d/ds(l_a^{ab}t_b)\mathrm{\Phi }`$ appearing in its normal variation (see Eq.(30)) will be non-vanishing, and it is no longer appropriate to discard a total derivative as we did in deriving Eq. (30). For consistency, we claim therefore that we need to modify Eq.(45) as follows:

$$𝐠𝐠l_a^{ab}t_b𝐧.$$
(64)
For Gauss-Bonnet, the second term reads $`2\beta K_{}𝐧`$. This mysterious term is precisely the tension associated with the edge energy $`_𝒞𝑑s\kappa `$. The projections of Eq.(45) along $`𝐥`$, $`𝐧`$ and $`𝐭`$ respectively then read

$`\alpha `$ $`[(KK_0)^22(KK_0)K_{}]2\beta K_{}^2`$ (65)
$`+`$ $`\sigma \kappa +\mu =0,`$
$`2`$ $`\alpha _{}K2\beta \dot{K}_{}+\sigma K_{}=0,`$ (66)
$`2`$ $`K_{}[\alpha (KK_0)+\beta K_{}]=0,`$ (67)
where we have used the fact that, at the edge,

$$\dot{𝐧}=K_{}𝐭+K_{}𝐥.$$
(68)
As was the case in the previous section, these coincide with the boundary conditions (60), (61), (62) when $`K_{}0`$.

## VI Conclusions

Whereas for a soap film, it is very simple to identify the forces operating on the edge, and so read off the boundary conditions on the bulk geometry, such an approach is less obvious for a membrane. However, we have demonstrated how simple geometrical and variational arguments may be exploited to derive the boundary conditions on the lipid membrane geometry. We have made no restrictive assumptions about the symmetry of the configuration. We then showed how these boundary conditions emerge from a balance of forces projected along a basis of vectors adapted to the edge.

We plan at a future date to examine the axially symmetric shape equation subject to the Gauss-Bonnet modified boundary conditions. Matching observed configurations with solutions of the shape equation in this more general setting provides a method, at least in principle, for determining the values of the two bulk rigidity moduli.

###### Acknowledgements.

RC is partially supported by CONACyT under grant 32187-E. JG and JS are partially supported by CONACyT under grant 32307-E and DGAPA at UNAM under grant IN119799.


text:
# Chapter 1 Bayesian Spectral Analysis of “MAD” Stars.

## Chapter 1 Bayesian Spectral Analysis of “MAD” Stars.

Nondas Sourlas11footnotemark: 1, David A. van Dyk 111Department of Statistics, Harvard University, One Oxford St., Cambridge, MA 02138, Vinay Kashyap22footnotemark: 2, Jeremy Drake22footnotemark: 2, and Deron Pease 222Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138

.

### 1 Overview

Getting reliable estimates of coronal metallicity ($`Z`$) from X-ray spectra obtained with instruments such as ASCA/SIS and Chandra/ACIS is very difficult, because the sole determinant of $`Z`$ is the ratio of line to continuum fluxes, which is not well-determined for low-resolution spectra. Here we propose new Bayesian methods which directly model the Poisson nature of the data. Our model also accounts for the Poisson nature of background contamination, blurring due to instrument response, and the absorption of photons in space. The resulting highly structured hierarchical model is fit using the Gibbs sampler, data augmentation, and Metropolis-Hasting. We demonstrate our methods with the X-ray spectral analysis of several apparently coronal metal abundance deficient (”MAD”) stars.

### 2 A Poisson Spectral Model

The model is designed to summarize the relative frequency of the energy of photons arriving at a detector. We model the photon counts in each bin as independent Poisson random variables. Specifically, we model the high energy tail of the ASCA spectrum (2.5-7.5 keV) as a combination of a Bremsstrahlung continuum and ten narrow emission lines, included at positions of known strong lines. This source model is combined with instrument response, the effective area of the instrument, and background contamination to model actual observed counts (van Dyk *et al.*, 2001).

Statistical analysis is based on two observations of the same source and two of the same background. We use sequential Bayesian analysis for the two source observations; the posterior distribution from the first analysis is used to construct an informative prior for the second. Non-informative priors were used in the first stage. Model fitting proceeds by using the EM algorithm to check for multimodality in the posterior. A Markov chain Monte Carlo algorithm is constructed to sample from that posterior. Three chains are used at each step of the analysis in order to assess convergence. Posterior inference is based on the second half of the draws of all the chains. The sensitivity of our results to the choice of prior was investigated by altering the prior. We used residual plots for both source and background to diagnose the fit of our models.

### 3 Results

We have measured the coronal metallicity for 4 stars, $`\alpha `$ Aur (Capella), $`\sigma `$ Gem, YY Gem, and II Peg. Our result for Capella ($`Z=0.73_{>0.24}^{<2.4}`$) are consistent with the independently determined value of $`Z=0.570.78`$ (Brickhouse *et al.*, 2000). Based on our analysis of $`\sigma `$ Gem ($`Z=0.81_{>0.59}^{<1.1}`$), we find that contrary to classical analyses ($`Z=0.25`$) it is not a MAD star (photospheric $`Z=0.6`$, Randich *et al.* (1994)). YY Gem is found to have sub-Solar metallicity ($`Z=0.46_{>0.23}^{<0.71}`$) but higher than EUVE/SW measurements ($`Z0.1`$; Kashyap *et al.* (1998)). For II Peg, we find $`Z=1.1_{>0.88}^{<1.3}`$, higher by a factor $`10`$ than preliminary results from Chandra/HETG; this discrepancy may be due to anomalies in the abundances of other elements, or to a high-temperature component to the emission measure.

Acknowledgments: The author gratefully acknowledges funding for this project partially provided by NSF grant DMS-01-04129 and by NASA contract NAS8-39073 (CXC). This chapter is a result of a joint effort of the members of the Astro-Statistics working group at Harvard University.