text:
# Strangeness Nucleation in Neutron Star Matter

## I Introduction

A ‘hot’ topic recently in nuclear physics has been the attempt to understand the production of strangeness in the expanding plasma thought to be created in relativistic collisions of heavy ions. There are several aspects to this complex and difficult problem, including: how (and, indeed, whether) a 3-flavor deconfined quark-gluon plasma is produced by the initial collision; how this strange matter equilibrates, expands, and cools; and how it eventually re-hadronizes, hopefully giving rise to a unique observable signal. The goal of these studies is to better understand the phase structure of QCD in the low-density, high-temperature region of the phase diagram.

A different but complementary problem involves the study of nuclear matter at high-density and low temperature, such as that existing in neutron star interiors. Here, matter consisting of neutrons, protons, and electrons (npe matter) is crushed during the gravitational collapse of the parent star. (Actually, what we call npe matter in this paper may also include muons and neutrinos.) As the star collapses weak interactions convert electrons and protons into neutrons, with electron neutrinos produced copiously as a by-product. Additionally, neutrinos of all flavors are created through brehmstrahlung-type collisions among the warm nucleons. As the neutrinos diffuse outward and radiate away, the proto-neutron star cools from an initial temperature of several tens of MeVs, and eventually settles into the familiar ground state of a neutron star with central density several times nuclear matter density.

It is now understood that at sufficiently high density, a phase transition will occur in which the neutrons and protons of npe matter are replaced by deconfined quarks (including strange quarks).  A related mechanism for the appearance of strangeness in neutron star matter is kaon condensation, which may become possible at densities somewhat lower than the density required for the deconfining transition described above. The possibility of kaon condensation was first pointed out by Kaplan and Nelson, who were motivated by the chiral Lagrangian prediction of an attractive interaction between $`K^{}`$ and nuclear matter.

Intuitively, one can understand the phenomenon of kaon condensation in the following way. As the density of npe matter is increased, and assuming that all particle species are in equilibrium with respect to the weak interactions, the electron chemical potential increases. Simultaneously, the effective mass of an in-medium $`K^{}`$ will decrease, due to the attrative interaction mentioned above. Therefore, at some critical density, the electron chemical potential (i.e., the energy of the electrons at the top of the Fermi sea) will become greater than the effective mass of the $`K^{}`$. It then becomes energetically favorable for the high-energy electrons to decay according to

$$e^{}\nu _e+K^{}$$
(1)
at which point a condensate of $`K^{}`$ will form.

This picture is somewhat over-simplified, as the production of kaons does not proceed exclusively by the decay of energetic electrons, but may also occur via fully-hadronic weak couplings, e.g., $`np+K^{}`$. Also, the notion of a uniform condensate of kaons is probably wrong. As pointed out by Glendenning , the presence of two separate conserved charges in neutron star matter (namely, baryon number and electric charge) gives rise to the possibility of a mixed phase of N-matter and K-matter existing over a wide range of densities and pressures in the star. (By “N-matter” we mean matter not containing kaons, but without the restriction of charge neutrality implied by “npe matter”.) This mixed phase will have the form of a Coulomb lattice – the now standard “pasta” structure consisting of droplets of the new phase immersed in a background of N-matter (at low densities), with rods and slabs replacing droplets as the density is increased, and finally with the role of the two phases interchanging in the transition toward homogenous K-matter as the density is increased still further.

In this paper, we will assume a model which predicts just this kind of first-order phase transition to kaon-condensed matter. For the description of nuclear matter, we use a Walecka-type model in which the interactions of strongly-interacting particles are mediated by $`\sigma `$, $`\omega `$, and $`\rho `$ mesons, treated in the mean-field approximation. Kaons are included in the model on the same footing as the nucleons, as described in the next section.

Our goal is, however, not to study the structure of the ground state predicted by this model (and its subsequent effect on the global neutron star properties such as the mass-radius relation), but rather to study the nucleation of kaons, that is, the process by which the initial proto-neutron star matter with zero strangeness acquires strangeness through the spontaneous appearance of K-matter droplets. We will begin this study by using the theory of homogenous nucleation due to Langer.  Here one assumes that at a given temperature, there are constant fluctuations producing small short-lived droplets of the new phase. Below a certain critical size, the surface energy cost of these droplets wins out over the volume term and the droplets shrink away. For some critical size, however, the free energy gained from the production of a large volume of the new (energetically favored) phase is just large enough to cancel the cost of surface energy, and the droplet will spontaneously grow. In the case studied here, because the two phases have non-zero electric charge densities, there is also a Coulomb term in the free energy, which becomes large for large droplets. Hence, a super-critical droplet will not grow forever, but will reach a stable size at which the energy gain between the volume and surface energies just balances the energy cost of the Coulomb energy.

Langer showed that the expected time needed for this kind of phase transition to proceed depended on the probability of a critical droplet being produced by thermal fluctuations, and also on the growth rate of such a critical droplet along the unstable direction in configuration space. Assuming a thermal distribution of these fluctuations, the nucleation rate per unit volume can be written

$$\mathrm{\Gamma }I_0e^{\mathrm{\Delta }F(R_{crit})/T}$$
(2)
where $`\mathrm{\Delta }F(R_{crit})`$ is the excess free energy of a critical droplet, $`T`$ is the temperature, and $`I_0`$ (the “prefactor”) is a microscopic fluctuation rate related to the growth rate of a super-critical droplet, the thermal conductivity of the medium, and other properties.  The prefactor is often approximated by dimensional analysis as simply $`I_0T^4`$. We will discuss this approximation later. Our goal is to apply this nucleation rate theory to the case of cooling proto-neutron star matter in order to better understand exactly when, where, and how the transition from npe-type matter to the kaon-condensed phase occurs.

Our work generally follows previous work on the nucleation of quark matter droplets in neutron star matter, which has been studied extensively.  The case of kaon condensation is physically unique, however, since it is the only realistic proposed example of a direct first order transition from npe-type neutron star matter to matter containing strangeness. In the case of the deconfinement transition mentioned above, for example, and even at densities, temperatures, and pressures for which 3-flavor quark matter is energetically favored, the transition is indirect in the following sense: an intermediate stage of 2-flavor quark matter will be produced first, with strange quarks then slowly appearing during a smooth crossover or second-order transition to the 3-flavor ground state. This latter transition is, of course, slow due to the weakness of the weak interactions. In this scenario, however, the original nucleation events which take one from npe matter to deconfined quark matter, need not involve the weak interactions at all.

The case of kaon condensation is radically different, since here there is no intermediate zero-strangeness state which might allow for fast nucleation followed by a slow but smooth growth of the strangeness-containing fields. Instead, the thermal fluctuations responsible for nucleation events must directly involve the weak-interaction processes which produce kaons. The difficulties posed by the widely varying timescales involved in these two sorts of processes (thermal fluctuations and weak interactions) will form the theme of our discussion.

The outline of the rest of the paper is as follows. In Section II we present the details of the nuclear mean field theory used for the subsequent discussion. In Section III we use this model to extract information about the free energy of K-matter droplets of different sizes, as the background baryon number density and temperature are varied. This allows us to use Eq. (2) to estimate the nucleation rate for the kaon-condensate phase transition. In Section IV we discuss in more detail the underlying mechanism for the fluctuations which produce the nucleation, and thereby analyze the trustworthiness of the estimates. Here we also make contact with the literature on Q-balls and argue that the problem of direct strangeness nucleation in neutron star matter may be profitably considered as an instance of Q-ball nucleation. Finally, in Section V we summarize the findings, discuss the relevance of our results to phenomenological issues in neutron star physics, and indicate some proposals for future investigation.

## II Mean-Field Theory Description of Kaon Condensation

In this section we briefly review the model proposed by Glendenning and Schaffner  which predicts a first-order transition from nuclear matter to the kaon-condensed phase. The model begins with a relativistic Walecka-type Lagrangian describing the neutron and proton fields, as well as the $`\sigma `$, $`\omega `$, and $`\rho `$ mesons which mediate their interactions:

$`_N=`$ $`\overline{\mathrm{\Psi }}_N(i\gamma ^\mu _\mu m_N^{}g_{\omega N}\gamma ^\mu V_\mu g_{rN}\gamma ^\mu \stackrel{}{\tau }_N\stackrel{}{R}_\mu )\mathrm{\Psi }_N`$ (5)
$`+{\displaystyle \frac{1}{2}}_\mu \sigma ^\mu \sigma {\displaystyle \frac{1}{2}}m_\sigma ^2\sigma ^2U(\sigma ){\displaystyle \frac{1}{4}}V_{\mu \nu }V^{\mu \nu }`$
$`+{\displaystyle \frac{1}{2}}m_\omega ^2V_\mu V^\mu {\displaystyle \frac{1}{4}}\stackrel{}{R}_{\mu \nu }\stackrel{}{R}^{\mu \nu }+{\displaystyle \frac{1}{2}}m_r^2\stackrel{}{R}_\mu \stackrel{}{R}^\mu ,`$
where $`m_N^{}=m_Ng_{\sigma N}\sigma `$ is the nucleon effective mass, which is reduced compared to the free nucleon mass due to the scalar field $`\sigma `$. The vector fields corresponding to the omega and rho mesons are given by $`V_{\mu \nu }=_\mu V_\nu _\nu V_\mu `$, and $`\stackrel{}{R}_{\mu \nu }=_\mu \stackrel{}{R}_\nu _\nu \stackrel{}{R}_\mu `$ respectively. $`\mathrm{\Psi }_N`$ is the nucleon field operator with $`\stackrel{}{\tau }_N`$ the nucleon isospin operator.

In addition to the usual kinetic, mass, and interaction terms for the nucleon and meson fields, the model also includes cubic and quartic self-interactions of the $`\sigma `$ field:

$$U(\sigma )=\frac{1}{3}bm_N(g_{\sigma N}\sigma )^3+\frac{1}{4}c(g_{\sigma N}\sigma )^4$$
(6)
where $`b`$ and $`c`$ are dimensionless coupling constants. These coupling constants (as well as the three nucleon-meson couplings: $`g_{\sigma N}`$, $`g_{\omega N}`$, and $`g_{rN}`$) are chosen to reproduce the empirical properties of nuclear matter at saturation density.

Kaons are included in the model in the same fashion as the nucleons, by coupling to the $`\sigma `$, $`\omega `$ and $`\rho `$ meson fields. There exist in the literature several meson-exchange Lagrangians which attempt to describe kaon-nucleon interactions. A detailed discussion of these models and their relation to the Chiral Lagrangian description proposed by Kaplan and Nelson  can be found in papers by Pons et al.  and Prakash et al. . In the present paper we employ the Lagrangian proposed by Glendenning and Schaffner.  The kaon Lagrangian is given by

$$_K=(𝒟_\mu K)^{}(𝒟^\mu K)m_K^2K^{}K$$
(7)
where $`K`$ denotes the isospin doublet kaon field. The covariant derivative $`𝒟_\mu =_\mu +ig_{\omega K}V_\mu +ig_{rK}\stackrel{}{\tau }_K\stackrel{}{R}_\mu `$ couples the kaon field to the vector mesons and the kaon effective mass term $`m_K^{}=m_Kg_{\sigma K}\sigma `$ describes its coupling to the scalar meson. $`\stackrel{}{\tau }_K`$ is the kaon isospin operator.

The vector coupling constants are determined by isospin and quark counting rules  and are given by $`g_{\omega K}=g_{\omega N}/3`$ and $`g_{rK}=g_{rN}`$. The scalar coupling is fixed by fitting to an empirically determined kaon optical potential in nuclear matter. The real part of this quantity has been determined from properties of kaonic atoms to lie in the somewhat wide range 80 MeV $`\stackrel{<}{}\text{ }U_K(n_o)\stackrel{<}{}\text{ }`$180 MeV . Here, we choose $`U_K(n_o)=120`$ MeV to fix $`g_{\sigma K}`$. Lower values of the kaon optical potential reduce the strength of the first-order transition and, eventually, produce instead a second-order phase transition, while higher values make the first-order transition stronger. While the issues discussed in this paper rely on the existence of a first-order transition, the qualitative conclusions are generally independent of the specific value of the coupling. So long as the transition is first-order, the main effect of changing $`U_K(n_o)`$ will be to change the critical density for the onset of the mixed phase, without severely affecting our discussion of the nucleation properties near this critical density.

The model as presented so far is a complicated, strongly interacting field theory which cannot be solved in any reasonable way. It is therefore standard to make a mean-field approximation, in which the meson field operators are replaced by their expectation values. Because of rotational invariance only the time-component of the vector fields $`V_\mu `$ and $`\stackrel{}{R}_\mu `$ can have a non-zero expectation value. Likewise, only the isospin 3-component of the isovector field $`\stackrel{}{R}_\mu `$ can be non-zero. The equations of motion for the meson (mean-) fields can be simply derived from the above Lagrangians and are given by:

$`m_\sigma ^2\sigma `$ $`=`$ $`{\displaystyle \frac{dU}{d\sigma }}+g_{\sigma B}(n_n^{(s)}+n_p^{(s)})+g_{\sigma k}m_K^{}f_K^2\theta ^2`$ (8)
$`m_\omega ^2\omega `$ $`=`$ $`g_{\omega N}(n_n+n_p)g_{\omega K}f_K^2\theta ^2(\mu _K+g_{\omega K}\omega +g_{rK}r)`$ (9)
$`m_r^2r`$ $`=`$ $`g_{rN}(n_pn_n)g_{rK}f_K^2\theta ^2(\mu _K+g_{\omega K}\omega +g_{rK}r)`$ (10)
where the meson fields $`\sigma ,\omega ,r`$ now represent the appropriate mean values. Here $`n_n`$ and $`n_p`$ represent the neutron and proton number densities, respectively, while $`n_n^{(s)}`$ and $`n_p^{(s)}`$ are the corresponding scalar densities. We have substituted $`K=(0,K^{})`$ and $`K^{}=\frac{1}{\sqrt{2}}f_K\theta e^{i\mu _Kt}`$ , where $`f_K`$ is the kaon decay constant and $`\theta `$ is a dimensionless kaon field strength parameter. We have neglected to write down the additional small contributions to the kaon-coupling terms due to finite temperature effects, though these are included in our code. (See Ref.  for a more explicit presentation of the details of this model for finite temperature.) We will be working in the bulk approximation (in which the meson fields do not vary with position) and so have set to zero the gradient terms which would otherwise appear in the equations of motion.

The equation of motion for the kaon field (in terms of $`\theta `$) is

$$0=\left((m_{K}^{}{}_{}{}^{})^2(\mu _{K}^{}{}_{}{}^{})^2\right)\theta $$
(11)
which indicates that the kaon effective mass $`m_K^{}=m_Kg_{\sigma K}\sigma `$ and the effective chemical potential $`\mu _K^{}=\mu _K+X`$ must be equal in order for the kaon field to take on a non-zero value. Here $`X=g_{\omega K}\omega +g_{rK}r`$ is the vector field contribution to the $`K^{}`$ energy. We also include a Langrangian describing non-interacting spin-$`\frac{1}{2}`$ particles to account for the presence of electrons, muons, and neutrinos.

The thermodynamic potential per unit volume for the nucleon sector is

$`{\displaystyle \frac{\mathrm{\Omega }_N}{V}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}m_\sigma ^2\sigma ^2+U(\sigma ){\displaystyle \frac{1}{2}}m_\omega ^2\omega ^2{\displaystyle \frac{1}{2}}m_r^2r^2`$ (13)
$`2T{\displaystyle \underset{i=n,p}{}}{\displaystyle \frac{\mathrm{d}^3k}{(2\pi )^3}\mathrm{l}n\left(1+e^{(E(k)\nu _i)/T}\right)}`$
where the nucleon energy $`E(k)=\sqrt{k^2+m_{N}^{}{}_{}{}^{2}}`$. The chemical potentials are given by $`\mu _p=\nu _p+g_{\omega N}\omega +\frac{1}{2}g_{rN}r`$ and $`\mu _n=\nu _n+g_{\omega N}\omega \frac{1}{2}g_{rN}r`$.

The other thermodynamic quantities can be calculated from $`\mathrm{\Omega }_N`$ in the standard way. The nucleonic contribution to the pressure, for example, is $`P_N=\mathrm{\Omega }_N/V`$, while the number densities and entropy densities are given by

$`n_{n,p}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_N/V}{\mu _{n,p}}}`$ (14)
$`=`$ $`2{\displaystyle \frac{\mathrm{d}^3k}{(2\pi )^3}f_F(E(k)\nu _{n,p})}`$ (15)
$`s_{n,p}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_N/V}{T}}`$ (17)
where $`f_F(ϵ)=\left(e^{ϵ/T}+1\right)^1`$ is the Fermi-Dirac distribution function. The neutron and proton scalar densities which enter in the equation of motion for the $`\sigma `$ field are given by

$$n_{n,p}^{(s)}=2\frac{d^3k}{(2\pi )^3}\frac{m_N^{}}{E(k)}f_F(E(k)\nu _{n,p})$$
(19)
The energy density is determined through the usual relation $`Ts_N=ϵ_N+P_N_i\mu _in_i`$ to be

$`ϵ_N`$ $`=`$ $`{\displaystyle \frac{1}{2}}m_\sigma ^2\sigma ^2+U(\sigma ){\displaystyle \frac{1}{2}}m_\omega ^2\omega ^2{\displaystyle \frac{1}{2}}m_r^2r^2`$ (21)
$`+2{\displaystyle \underset{i=n,p}{}}{\displaystyle \frac{\mathrm{d}^3k}{(2\pi )^3}E(k)f_F(E(k)\nu _i)}`$
$`+`$ $`{\displaystyle \underset{i=n,p}{}}n_i(\mu _i\nu _i).`$ (22)
The thermodynamic potential for the lepton species present (electrons, muons, and neutrinos) is given by

$$\frac{\mathrm{\Omega }_L}{V}=\underset{l}{}Tg_l\frac{\mathrm{d}^3k}{(2\pi )^3}\left[\mathrm{l}n\left(1+e^{(E_l(k)\mu _l)/T}\right)+\mathrm{l}n\left(1+e^{(E_l(k)+\mu _l)/T}\right)\right]$$
(23)
where $`\mu _l`$ denotes the chemical potential for lepton species $`l`$, and $`g_l`$ are the spin-degeneracies: $`g=2`$ for electrons and muons, $`g=1`$ for neutrinos. $`\beta `$-equilibrium requires the following constraints on the chemical potentials:

$$\mu _K=\mu _e\mu _{\nu _e}=\mu _n\mu _p=\mu _\mu \mu _{\nu _\mu }$$
(24)
The lepton contributions to the pressure, energy density, and entropy are determined from Eq.(23) in the usual way.

Finally, the thermodynamic potential for the kaons is given by

$`{\displaystyle \frac{\mathrm{\Omega }_K}{V}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}f_K^2\theta ^2\left[(m_K^{})^2(\mu _K+X)^2\right]`$ (26)
$`+T{\displaystyle \frac{d^3p}{(2\pi )^3}\mathrm{l}n\left(1e^{(\omega ^{}(p)\mu _K)/T}\right)}`$
where $`\omega ^{}(p)=\sqrt{p^2+(m_K^{})^2}X`$ is the in-medium energy of a $`K^{}`$ with momentum $`p`$. Again, the kaon contribution to the other thermodynamic quantities can be found by differentiation. For a detailed derivation and presentation of the thermodynamics of this model, see Ref. . For discussions of similar studies involving mixed-phases in mean-field-theories in several different physical contexts see, e.g., .

Our goal here is to study the rate at which K-matter nucleates in a background of npe-matter at various temperatures and densities. As discussed previously, this involves calculating the free energy of droplets of various sizes, in particular, the free energy of the critical droplet configuration. The first step toward this end is to produce a description of a droplet of K-matter of arbitrary radius. This involves solving the meson field equations of motion subject to various constraints.

Generally, one solves simultaneously two versions of the meson field equations above: one with $`\theta =0`$ describing the N-matter component of the mixed phase, the other with $`\theta 0`$ (with the specific value of $`\theta `$ determined self-consistently through Eq. (11)) describing the K-matter. One requires chemical equilibrium between the two phases, i.e., that the relevant chemical potentials in the two phases match. Additionally, the requirement of overall electric charge neutrality means that the electric charge densities in the two phases must be of opposite sign. One then uses the respective charge densities to calculate $`\alpha =q^{(N)}/(q^{(K)}q^{(N)})`$, the volume fraction of the K-matter phase. (Here $`q^{(N)}`$ and $`q^{(K)}`$ are the charge densities of the K-matter and N-matter, respectively.) Because we are concerned with droplets of small radius, it is also crucial to include correct mechanical equilibrium between the two phases. In Ref. this requirement was shown to affect rather dramatically the bulk properties of the two phases. (See also .) For a spherical droplet of K-matter in a background of N-matter, this constraint reads:

$$P_KP_N=2\sigma /R$$
(27)
where $`\sigma `$ is the surface tension between the two phases, and $`R`$ is the radius of the droplet.

Finally, it should be noted that we work at constant baryon number density. That is, in constructing a sequence of droplets of varying radii, we require that the overall baryon density $`n_B=\alpha n_B^{(K)}+(1\alpha )n_B^{(N)}`$ be held fixed at some specified value. Here $`n_B^{(N,K)}`$ are the baryon number densities of the N- and K-phases.

Once the equations of motion are solved self-consistently subject to these constraints, it is possible to calculate the free energy of a given droplet. The total bulk energy density of a given droplet configuration is defined analogously to the overall baryon number density:

$$ϵ_{bulk}=\alpha ϵ^{(K)}+(1\alpha )ϵ^{(N)}$$
(28)
where the energy density in each phase includes the contributions from the nucleons, leptons, and kaons, as appropriate. The total bulk energy is then found by multiplying this energy density by the total volume of a Wigner-Seitz cell (defined as the spherical region including one droplet and containing zero total electric charge), $`V_{WS}`$. This energy is then supplemented by the surface and Coulomb contributions:

$$E(R)=ϵ_{bulk}V_{WS}(R)+4\pi R^2\sigma +E_{Coul}$$
(29)
where $`E_{Coul}`$ is found by integrating the electric field energy density $`ϵ_{Coul}(r)=\frac{1}{2}|E(r)|^2`$ throughout the Wigner-Seitz cell. (Here $`|E(r)|`$ is the magnitude of the electric field determined by Gauss’ Law.) The total entropy of a droplet configuration is found analogously, by volume-fraction-weighted averaging of the individual phase entropy densities, then multiplying by the volume of the Wigner-Seitz cell. Thus, the total free energy of a droplet can be calculated as:

$$F(R)=E(R)TS(R).$$
(30)
Only the R-dependence has been indicated explicitly, but, of course, the energy and entropy both depend strongly on the fixed baryon number density, the temperature, neutrino fraction, etc.

We are interested, however, not merely in the free energy of a given configuration, but, rather, the change in free energy required to produce various droplets. Hence, we also solve for pure, charge neutral npe matter at a given baryon density and calculate its free energy density, $`f_{npe}=ϵTs`$. We are then led to define

$$\mathrm{\Delta }F(R)=F(R)f_{npe}V_{WS}.$$
(31)
which represents the free energy cost of a transition from electrically neutral npe matter (the initial state of an evolving proto-neutron star) to a single droplet of kaon-condensed matter of radius $`R`$. Fig. 1 shows $`\mathrm{\Delta }F(R)`$ as a function of $`R`$ near the critical density for the transition, for two different values of the surface tension, $`\sigma =20MeV/fm^2`$ and $`\sigma =30MeV/fm^2`$.

We have shown the free energy curves for two values of $`\sigma `$ in order to illustrate the role this quantity plays. As expected, smaller values of $`\sigma `$ reduce both the size and free energy cost of a critical droplet. In what follows, we will simply pick the value $`\sigma =30\mathrm{M}eV/fm^2`$, a value suggested by Glendenning’s study of the boundary between the two phases , as well as by our own earlier work on this model .

In the following section we will use this type of curve to acquire information about the free energy cost of a critical droplet under various conditions of temperature, density, and neutrino fraction.

## III Droplet Free Energy and Nucleation Rates

As mentioned in the introduction, one can understand the production of kaons in the neutron star matter as being due to the weak decay of an electron (or, equivalently, the change of a neutron into a proton plus a kaon). By definition, above the critical density this transition is energetically favorable, i.e., exothermic. So we may schematically write

$$e^{}\nu _e+K^{}+\mathrm{H}eat.$$
(32)
The excess heat generated by the production of kaons will eventually be radiated away in the form of photons and neutrinos as the neutron star cools, but that is not our interest here. Rather, our goal is to understand how the free energy of droplet configurations depends on the temperature. Qualitatively, one can guess the correct answer by applying LeChatelier’s principle to the equilibrium indicated in Eq.(32): raising the temperature will tend to push the equilibrium back toward the left. That is, raising the temperature should lower the energetic favorability of kaon droplets.

This prediction is borne out by our calculations, as illustrated in Fig. 2. As the temperature is increased from zero, the free energy of kaon droplet configurations increases relative to that of neutral npe matter at the same density and temperature. This increase (which is as large as several hundred MeV for typical droplets) means that, at a density where the mixed phase would be the ground state at $`T=0`$, the mixed phase is no longer favored at higher temperature. In other words, turning up the temperature increases (slightly) the critical density for the transition.

However, because the kaon phase droplets must be produced by thermal fluctuations, we expect the thermal nucleation rate of the kaon phase to be much larger at higher temperatures. In order to estimate the speed with which droplets of the kaon phase are produced, we use the Langer nucleation rate theory discussed in the introduction. The nucleation rate per unit volume is estimated as:

$$\mathrm{\Gamma }T^4e^{\mathrm{\Delta }F(R_{crit})/T}$$
(33)
where T is the temperature and $`R_{crit}`$ is the critical droplet radius (at a given temperature and baryon density). The free energy cost of this critical droplet configuration (i.e., the height of the barrier which must be crossed to produce a stable droplet of the new phase) can be easily read off of graphs like the ones already shown. One may then plug directly into Eq. (33) to give the nucleation rate per unit volume. In order to convert this rate into an intuitively meaningful quantity, we calculate the expected time for a single nucleation event in a single typical Wigner-Seitz cell of volume $`V_{WS}=10^3fm^3`$. This time is given by

$$\tau =\frac{1}{\mathrm{\Gamma }V_{WS}}=\frac{e^{\mathrm{\Delta }F(R_{crit})/T}}{V_{WS}T^4}$$
(34)
In Fig. 3 we show this nucleation time as a function of density for several different temperatures. (Note the log scale!) As expected, the nucleation time is a very strong function of the temperature. At a temperature of $`0.1`$ MeV, the expected nucleation time is many, many times longer than the age of the universe across the entire density range of the mixed phase. At $`T=1.0`$ MeV the nucleation time is prohibitively long at the lower end of the mixed phase density regime, but is less than one second at densities above about $`n_B0.55fm^3`$, which is something like $`\frac{1}{2}n_o`$ above the critical density. At higher temperatures, the nucleation proceeds almost immediately across the entire density range.

Naively, this leads to the conclusion that the neutron star settles into its ground state without any delay as it cools down from an initial temperature of several tens of MeV. In the initial, hot conditions, thermal fluctuations are sufficiently fast and sufficiently numerous to seed K-matter droplets wherever those droplets are energetically favored. The sizes and distances between adjacent kaon structures may undergo tiny changes as the matter cools, and perhaps the outer edge of the mixed phase extends outward somewhat as the critical density decreases, but generally the neutron star includes the full mixed phase from birth.

This picture is complicated by at least two factors which have not been discussed explicitly until now. The first of these is the presence, in the early stages of proto-neutron star cooling, of a significant neutrino fraction. The second question is whether or not we should believe the rate estimates just given, since the nucleation rate prefactor (estimated above as $`T^4`$) describes the rate of microscopic fluctuations, which, in the present case, consist of weak interaction processes. These processes are, after all, weak, so one might doubt that the naive estimate based simply on dimensional analysis is appropriate. This issue will be discussed in a subsequent section; for now, we will turn to the question of the effects of neutrinos on the nucleation of K-matter. (For a related discussion see Ref.)

Returning to Eq.(32) and again applying LeChetalier’s principle, we guess that the presence of a non-zero density of electron neutrinos will (like high temperature) suppress the appearance of K-matter. This guess turns out to be correct, but for slightly complicated reasons. Turning on a non-zero electron neutrino fraction

$$Y_{\nu _e}=\frac{n_{\nu _e}}{n_B}$$
(35)
forces the chemical potentials for electrons, kaons, and the baryons to adjust according to the constraint of constant $`n_B`$ and $`\beta `$-equilibrium, Eq. (24). Surprisingly, doing this actually lowers the free energy cost of kaon droplets, as indicated in Fig. 4.

The reason for this can be understood as follows. The main effect of turning on a chemical potential for neutrinos, is to increase as well the chemical potential for electrons since the constraint of fixed baryon number density doesn’t allow $`\mu _n`$ and $`\mu _p`$ much freedom to adjust. But increasing the number density of electrons decreases the (positive) electric charge density of the N-matter outside a kaon droplet. Hence, more of this matter is needed to cancel the negative charge of the droplet itself, and the volume fraction of K-matter, $`\alpha `$, decreases. But the overall free energy density depends on $`\alpha `$ through Eq.(28). Finally, since the energy density of K-matter is somewhat larger than that for electrically neutral npe matter (due not only to the presence of kaons, but also to the higher local baryon density) while that of the surrounding N-matter is somewhat lower than neutral npe matter (due to the correspondingly lower baryon number density here), the overall energy density of the mixed phase is actually decreased by the decrease in $`\alpha `$ coming from the non-zero neutrino density.

However, because of the requirement of overall electric charge neutrality and the fact that the K-matter is always negatively charged, the N-matter component outside of the droplet must have a positive electric charge density. As we have just seen, the presence of even a small non-zero chemical potential for electron neutrinos also increases the density of electrons, and thereby reduces the positive charge in this region. As expected, a larger neutrino fraction will decrease the electric charge density here even more, and eventually a point is reached at which the charge density ceases to be positive. Then it is no longer possible to define a Wigner-Seitz cell, i.e., no longer possible to satisfy the requirement of electric charge neutrality.

Actually, this point is reached for very moderate values of the neutrino fraction, especially at low densities (near the critical density for the transition) where, by definition, the electric charge density in the N-matter phase is positive, but small. We find that the pressure equilibrium condition Eq.(27) also has a comparable effect on the electric charge densities via the electron chemical potential. (This issue was discussed in Ref..) In particular, larger structures at a given baryon density have lower electron chemical potentials; hence it is possible to elude the effect of non-zero $`Y_\nu `$ by producing larger K-matter droplets. Thus, at a given $`Y_\nu `$ and a given baryon density, there will be a minimum size $`R_{min}`$ consistent with charge neutrality, and such structures are therefore not allowed. Electron Debye screening will also play a role in forbidding structures with large radii, so the region of $`n_BY_\nu `$ parameter space in which no mixed phase can exist is even somewhat larger than suggested here.

We show in Fig.5 this minimum radius as a function of baryon density, for several different values of the electron neutrino fraction, $`Y_\nu `$. At a given $`Y_\nu `$, the minimum radius allowed by global charge neutrality begins to diverge as one comes down in density. The density at which this quantity diverges (or, if we were to consider also the effects of Debye screening, the density at which this quantity exceeds 5-10 fm) acts as an effective critical density for the onset of the mixed phase. Hence, early in the evolution of the proto-neutron star (PNS), when the neutrino fraction is (at least) a few percent, no mixed phase will be produced at densities less than $`0.600.65fm^3`$. This is to be compared to the critical density at zero temperature and zero neutrino fraction, $`0.49fm^3`$. Actually, the electron neutrino fraction in the first seconds of PNS evolution may be closer to ten percent, in which case the effective critical density below which no stable kaon mixed phases exists will be pushed upward to around $`0.80fm^3`$, several factors of $`n_0`$ above the nominal ($`T=0`$, $`Y_\nu =0`$) critical density. The upshot is that, in the neutrino-rich conditions of the early PNS evolution, the critical density for the onset of the kaon mixed phase is increased significantly compared to the nominal critical density.

This obviously modifies the naive inference from Fig.3 that the mixed phase will be produced quickly and easily during the early, hot part of the PNS evolution. Because the high temperatures occur at a time of high neutrino density, there is a wide range of densities over which no mixed phase can be formed at these early times. Hence, due to the suppression of K-matter by neutrinos, a PNS with central density not too far above the nominal critical density for the onset of the mixed phase may fail to produce any K-matter during its cooling, until the temperature is so low that nucleation can not occur in reasonable times. This suggests a scenario in which the neutron star could exist indefinitely in a metastable state. However, further analysis is required to support this hypothesis – in particular, we must address the question of the reliability of the nucleation rates in Fig.3 given the problem of simultaneous weak interactions.

## IV Fluctuations and the Problem of Simultaneous Weak Interactions

As a first guess, one might suppose that the nucleation of K-matter droplets is driven by thermal fluctuations in the local number density of baryons. At a sufficiently high baryon number density, the phase transition to the kaon condensed phase becomes second order, and the kaon field may grow smoothly without having to overcome an energy barrier. Hence, one might suppose that the first-order transition could occur by this same mechanism operating locally: a thermal fluctuation produces a baryon overdensity in a local region of radius $`12fm`$. The kaon field then spontaneously “fills in” this region, producing a stable droplet of K-matter which could then grow to the stable size. In this picture, the nucleation rate would be governed by the frequency of sufficiently over-dense and sufficiently large baryon number fluctuations.

However, there is an immediate and fatal problem with this basic picture. Sufficient fluctuations in baryon density are no-doubt plentiful. A typical $`12fm`$ region contains of order ten baryons, so assuming simple $`\sqrt{N}`$ fluctuations, a factor of two overdensity is only three standard deviations away from the mean. But the lifetime of such a fluctuation is limited to strong interaction timescales. Indeed, in nuclear matter at these densities, the speed of sound is a sizable fraction of the speed of light, $`c`$. Hence, the lifetime of a baryon number density fluctuation is of order $`\tau R/c10^{23}sec`$. As is well-known from studies of a second-order kaon condensate phase transition, however, the time-scales needed for the development of an appreciable value of the kaon field are $`1015`$ orders of magnitude longer than this.  Similar conclusions are found for the appearance of strange quark matter from an initial state of 2-flavor quark matter in neutron stars.  That is, not surprisingly, the time-scale for the development of the kaon field is typical of the weak interactions.

The production of K-matter droplets, therefore, cannot be driven by density fluctuations in the background of baryons. The microscopic kaon-production rate is simply too slow to keep up with the strong interaction timescales that govern such fluctuations, so long before any appreciable value for the kaon field in the fluctuation has developed, the fluctuation evaporates.

We conclude that the nucleation rate will be goverened by thermal fluctuations producing fluctuations in the local kaon number density itself. Moreover, so long as the timescale for the spontaneous growth of a critical droplet is small compared to weak interaction timescales, this rate can depend only on the thermal number density of kaons,

$$n_K^{(thermal)}=\frac{d^3p}{(2\pi )^3}f_B(\omega ^{}(p)\mu _K)$$
(36)
where $`f_B`$ is the Bose-Einstein occupation probability and $`\omega ^{}(p)`$ was defined in Sect.II. We have shown in Fig.6 this thermal kaon number density for several values of the temperature and neutrino fraction in npe-type neutron star matter as a function of density.

It is then simple to estimate the probability of a sufficient “seed” of kaons appearing. This calculation is reminiscent of the standard undergraduate statistical mechanics problem of calculating the probability that all of the atoms in a box of gas are found to be in a certain small region of the box. Here, we must calculate the probability that a certain number $`N`$ of the kaons in a kaon gas of density $`n_K^{(thermal)}`$ spontaneously appear in some small region of space with volume $`V`$.

For definiteness, consider a large box of volume $`V_0`$ and total kaon number $`N_0`$, with $`N_0/V_0=n_K^{(thermal)}`$. If we treat the kaons as classical particles, then the probability of $`N`$ kaons being found in a small sub-region $`V`$ of $`V_0`$ is given by the binomial theorem:

$$P(N)=\frac{N_0!}{N!(N_0N)!}\left(\frac{V}{V_0}\right)^N\left(1\frac{V}{V_0}\right)^{(N_0N)}$$
(37)
Assuming $`VV_0`$ and $`NN_0`$, we may approximate $`N_0!(N_0N)!N_0^N`$ and neglect $`N`$ in the exponent $`(N_0N)`$. This gives

$$P(N)\frac{1}{N!}\left(\frac{N_0V}{V_0}\right)^N\left(1\frac{V}{V_0}\right)^{N_0}.$$
(38)
But $`N_0V/V_0=<N>`$, where $`<N>=Vn_K^{(thermal)}`$ is the average number of thermal kaons in the volume $`V`$. Using $`<N>/V=N_0/V_0`$ we have finally

$`P(N)`$ $``$ $`{\displaystyle \frac{<N>^N}{N!}}\left(1{\displaystyle \frac{<N>}{N_0}}\right)^{N_0}`$ (39)
$``$ $`{\displaystyle \frac{<N>^N}{N!}}e^{<N>}.`$ (40)
(In practice, $`n_K^{(thermal)}1fm^3`$ so the exponential in the last line above can be ignored.)

The number of kaons in a critical droplet varies somewhat with density. Near the critical density (at zero temperature and zero neutrino fraction) the radius of a critical droplet is a few $`fm`$, for a critical kaon number of order 100. At higher densities, the critical radius drops to only $`12fm`$, and the corresponding critical kaon number drops to only a few. In Fig. 7 we plot the time needed to form such a critical seed (for various $`N`$) as a function of thermal kaon number density. To estimate the seed-production time from the probability considered above, we multiply each probability by the maximum possible number of discrete seeding locations within a typical Wigner-Seitz cell. This is of order $`V_{WS}n_B1000fm^3.5/fm^3500`$. We then assume that the arrangement of thermal kaons is “reshuffled” on the fastest possible time-scale, say $`\tau 10^{23}sec`$. (Note that this assumption is extremely generous. The actual reshuffling time must be at least several orders of magnitude greater than this, but here a generous upper limit on the seed-production rate will suffice.) Then the typical time needed for a seed of $`N`$ kaons to be produced by this method is given by

$$\tau (N)\frac{10^{23}sec}{500P(N)}$$
(41)
with P(N) given as a function of $`n_K^{(thermal)}`$ above. Actually, we should replace $`P(N)`$ here with $`_{nN}P(n)`$ but in practice each term in the sum is negligible compared to the one previous, so the difference will not affect the result.

We see that the number $`N_{critical}`$ is indeed critical to the determination of the rate. If ten (or more) kaons are needed to produce a critical droplet, the probability of a sufficient number all showing up in the same place at the same time is extremely small, and one must wait prohibitively long for a seed to ever be produced, especially at lower temperatures where $`n_K^{(thermal)}`$ is extremely small. At higher densities where $`N_{crit}`$ is lower, of order 3-5, the probability is much higher the time needed for a seed to be produced by random reshuffling may be of order seconds or less for reasonable values of the temperature.

Comparing Fig. 6 and Fig. 7, we see that the time needed for a seeding event during the initial hot proto-neutron star configuration ($`T=30MeV`$, $`Y_\nu =0.03`$) is on the order of a few seconds or less if $`N_{crit}510`$. However, as we see from Fig. 5, electric charge neutrality forbids droplets smaller than $`2fm`$ for densities below $`0.65fm^3`$. This is well above the ($`T=0`$, $`Y_\nu =0`$) critical density, and it is likely that such a high density is not reached in the proto-neutron star core, especially during the early evolution when the star is still hot and relatively rarefied.

After a few tens of seconds, the PNS temperature drops to $`T=10`$ MeV with the neutrino fraction also dropping to, say, $`Y_\nu =0.01`$. Looking again at Fig. 6, we see that the number density of thermal kaons has dropped dramatically due to the lower temperature. The electric charge neutrality constraint is relaxed somewhat, with droplets of radius $`12fm`$ now being allowed at densities below about $`n_B0.58fm^3`$. There is a possibility, considering Fig. 7, that seeding may occur with reasonable speed for densities above this value, where $`N_{crit}`$ is of order $`35`$. But again, there is a fairly rigid effective critical density below which nucleation cannot occur. As the star continues to cool, the number density of thermal kaons drops dramatically, and the seeding time becomes prohibitively long at all but the very highest densities, $`57`$ times nuclear matter density, reasonably thought to exist in the neutron star core.

The crude estimates made here for seed production in a homogenous but fluctuating background of thermal kaons are similar to the estimates in  of the growth process for sub-critical Q-balls in a scalar field with a conserved charge. These authors assume a thermal distribution of sub-critical Q-balls with relatively large N, and then consider the rate at which additional charge is accreted to (and released by) the Q-ball via a random-walk process in order to estimate the rate at which critical droplets are produced. In our case, the number of kaons in a critical droplet is not much greater than unity, so we use the Poisson statistics of Eq. (40) rather than a Gaussian distribution with $`\sqrt{N}`$ fluctuations.

Actually, there is an exact analogy between the physics of kaon nucleation and the problem of Q-ball nucleation, i.e., the decay of metastable field configurations in the presence of a conserved charge.  The effective theory for the kaon field in our model will be precisely the theory of a complex scalar with an effective potential supporting non-topological solitons, i.e., stable droplets of K-matter. In the limit of infinitely slow weak interactions, this effective theory will contain a global $`U(1)`$ flavor symmetry corresponding to the conservation of strangeness. In this limit, the seeding of K-matter droplets will occur via a process formally identical to that described in Ref. , namely, a small uniform initial charge density (in our case, the background of thermal kaons) producing the required “bounce” configuration by flowing toward a seeding center. As discussed in , the action associated with this configuration will be significantly greater than the corresponding case for a real scalar which carries no conserved charge. The typical nucleation rates in the case of a charged scalar will thus be exponentially slower, due to the need for a global redistribution of charge.

The two estimates of nucleation rates in the present paper can be thought of as limits in the cases of infinitely slow and infinitely fast weak interactions. In reality, of course, the weak interactions proceed at some intermediate speed. In the effective theory of the kaon field described above, inclusion of strangeness-changing weak reactions will correspond to the addition of a small symmetry-breaking term in the effective potential, e.g.,

$$U_{eff}(|K|)U_{eff}(|K|)+\lambda K$$
(42)
where $`\lambda `$ is a strongly temperature-dependent constant that characterizes the rate at which kaons can be produced/absorbed by background scattering processes involving nucleons and leptons.

The actual rate of K-matter nucleation will depend on both types of symmetry breaking, that is (1) the initial background thermal kaon density and (2) the dynamical production of new kaons, though it is not obvious which mechanism will dominate the nucleation for physically relevant temperatures. In order to treat this problem reliably, one must construct realistic bounce configurations in the kaon effective theory as the initial background density and dynamical breaking parameter $`\lambda `$ are varied with temperature, impose imaginary-time periodicity corresponding to the (inverse) temperature, and calculate the action of the bounce. Qualitatively, one expects the resulting nucleation times to be intermediate between the fast-weak-interactions and no-weak-interactions limits in the current paper. At very least, one should expect that the slowness of the weak interactions should modify the nucleation rate prefactor $`I_0`$ to a more realistic dimensional analysis estimate

$$I_0G_F^2T^8\frac{T^8}{M_W^4}$$
(43)
where $`M_W`$ is the mass of the weak gauge bosons. At a temperature of $`10`$ MeV, this estimate would reduce the expected nucleation time by a factor of approximately $`10^{16}`$ from the times shown in Fig. 3. As mentioned above, the issue of nearly-conserved strangeness will also affect the exponential factor in the rate equation, thus suppressing the rate even further and perhaps bringing the results closer to the estimates in Fig. 7. Only the reliable calculation outlined above (and currently underway by the present author) will answer this issue with any certainty, however.

## V Discussion and Conclusions

We have, then, an intriguing possible scenario in which the PNS manages to settle into its ground state without forming the kaon-condensate mixed phase which is the true ground state of the system. At high initial temperatures, the nucleation of K-matter droplets is suppressed by the presence of neutrinos, even though at these high temperatures the seeding of droplets is in principle fast. Over a wide range of densities kaon droplets are produced copiously by thermal fluctuations, but they are not yet stable due to the increase in effective critical density. As the star cools, the restriction coming from the presence of neutrinos is relaxed, but the intrinsic fluctuation rate drops.

For an initial PNS core density that is not too far above the nominal ($`T=0`$) critical density for the formation of a K-matter mixed phase, therefore, it is likely that the star will cool not into the true ground state, but, rather, into a meta-stable state consisting of electrically neutral npe matter. This scenario may be relevant to understanding various phenomenological issues.

For example, the apparent existence of anomalously heavy neutron stars with masses $`M2M_{sun}`$  might be explained by the anomalously stiff equation of state of npe-type matter relative to matter with a kaon condensate. Generally, if the various possible phase transitions thought to occur in neutron star matter can be avoided by the impossibility (or, equivalently, extreme slowness) of nucleation of the new phase, a relatively stiff equation of state may sometimes be maintained over a more extended range of density than would be expected naively. Hence, the existence of such heavy stars may not be sufficient evidence to rule out the existence of kaon-condensation at $`35n_0`$, especially if these stars were born with smaller masses and only subsequently (that is, once cold) acquired larger masses via accretion from companion stars.

Additionally, metastability of the sort introduced above may potentially be useful in understanding the properties of GRB’s or other poorly-understood explosive events. A relatively light PNS, as we have argued, may cool and deleptonize without the K-matter mixed phase forming, even when the star’s core density exceeds the critical density for the transition. As is evident from Fig. 6, however, the thermal kaon density increases monotonically and steeply with density, so that, even at the very low temperatures $`T1`$ MeV eventually attained in the neutron star, there is some density at which seeding may become possible in a reasonable time. At very least, with increasing density, one eventually encounters the second-order point at which the kaon field may be produced smoothly with no need for seeding. Hence, if an initially metastable neutron star begins to accrete matter from a companion binary (or, additionally, if an initially rotating metastable neutron star gets spun down via accretion) the central density may increase sufficiently for K-matter to begin to appear.

Once the kaon matter appears, however, there will be a feedback effect, due to the softening of the equation of state. The production of a small quantity of K-matter in the core would allow the star to contract slightly, thus increasing the density in the core, and increasing the size of the region in which kaonic matter can appear. Further kaon production leads to further collapse, and vice versa. Thus if mass accretion and/or spin down results in the critical density for the onset of a second-order kaon-condensate transition being reached in the neutron star core (or, equivalently, if one reaches the low-temperature effective critical density at which a mixed phase can be nucleated spontaneously with sufficient speed) one would expect an explosive event in which the star contracts significantly, resulting in the release of a tremendous amount of energy (of order tenths of $`M_{sun}`$). Cheng and Dai have discussed a similar proposal in which accretion-induced conversion to strange quark matter is suggested as a possible explanation for Gamma-Ray Bursters.

One especially interesting aspect of such a collapse is its potentially turbulent nature. The picture is of a pure K-matter core seeding the mixed phase through several kilometers of material above it. In effect, the K-matter boils off of the outer edge of the second-order core and floats upward to form the mixed phase throughout the entire region of the mixed phase’s energetic favorability. This implies an upward and downward transfer of matter that closely resembles turbulent convection, but in which strangeness rather than heat is the substance being convected.

As mentioned at the end of the previous section, more reliable calculations need to be performed in order to better understand how the slowness of the weak interactions affect the original nucleation rate estimates based on Langer’s formula. We have argued that the correct framework for these future calculations involves the formalism of quantum tunneling (or thermal activation) in a theory with a (nearly) conserved global charge representing strangeness. It is also worth mentioning that the scenario outlined here is the most realistic potential application of the Q-ball nucleation formalism developed in ; theorists up to now have relied on supersymmetric models to find possible theories containing charged scalars supporting non-topological solitons.

Acknowledgements

Sanjay Reddy, Eduardo Fraga, and Guy Moore are acknowledged for helpful discussions, though any inaccuracies in the present paper are purely the responsibility of the author. This work is supported in part by a National Science Foundation Graduate Research Fellowship and by the US Department of Energy grant DE-FG03-00ER41132.


format: text

text:
# Does the melting behavior in ultrathin metallic nanowires begin from the surface ?

## Abstract

Surface melting temperature is well-known significantly lower than the bulk melting point. But we find that the interior melting temperature in ultrathin nanowires is lower than that of the surface melting . The thermal stability of helical multi-walled cylindrical gold nanowires is studied using molecular dynamics simulations. The melting temperature of gold nanowires is lower than the bulk value, but higher than that of gold nanoclusters. An interesting interior melting is revealed in the gold nanowires and the thermodynamics is closely related to interior structures. The melting starts from the interior atoms, while the surface melting occurs at relatively higher temperature. We propose that the surface melting represents the overall melting in ultrathin metallic nanowires.

The melting behavior of nanoparticles and nanorods have been demonstrated dramatically different from the bulk both experimentally and theoretically . The melting process of a crystalline starts from the surface layer and propagates into the interior. Thus, the surface melting temperature is significantly lower than the bulk melting point. Similarly, people may ask whether the surface melting temperature is lower than the overall melting temperature in clusters and nanowires. Berry considers that “dynamic coexistence” or surface melting happens in the melting process of small clusters before the total melting. For crystalline nanowires, Tosatti found that surface melting temperature for Pb wires was also lower than the total melting temperature. Experimentally, surface melting is involved in the melting process of nanoparticles and nanorods . Especially, Schmidt found a broad peak for heat capacity of Na$`{}_{}{}^{+}{}_{139}{}^{}`$ clusters, maybe implying the coexistence of solid-like and liquid-like phases before the total melting . Two major effects are responsible for these different melting behavior in nanoparticles and nanorods. One is the large surface-to-volume ratio in these nanostructures. The other is quantum confinement effect in finite size systems. Surface atoms have fewer nearest neighbors and weaker binding, which may induce an earlier surface melting behavior. On the other hand, a close relationship between the melting and the structure features is found in the cluster.

Recent studies demonstrated that ultrathin metallic nanowires have quite different structural and properties from those of bulk, clusters and crystalline nanowires . The helical multi-walled cylindrical structures have been widely found in 1-3nm size range of metallic nanowires both experimentally and theoretically . This kind of novel structure will bring some bizarre melting features different from the above-mentioned systems. To our knowledge, less efforts are focused on their thermodynamics so far, although this kind of ultrathin metallic nanowires have attracted great interests .

Furthermore, the ultrathin nanowire has some characteristics similar to the cluster, crystalline wire and bulk. It may provide an opportunity for comprehensively understanding these types of matter states and their relations. For example, the surface-to-volume ratio is a non-zero value in these nanostructures, while it approaches toward zero in the bulk. Surface and core atoms are expected to play different role during the melting process. However, it is difficult to definitely distinguish the overall melting and surface melting of clusters experimentally because the cluster’s signal spectrum will not exist when the surface melting occurs. Therefore, we can employ the ultrathin gold nanowire as a representation to explore these problems.

In this letter, the thermal stabilities of gold nanowires with helical multi-walled cylindrical structures are studied by using molecular dynamical (MD) simulations. We start from the optimized structures from previous works , which were imaged by electron microscopy. The interaction between gold atoms is described by a glue potential  and the periodic boundary condition is applied along the wire axis direction to model the wire with sufficient length wire. The length of supercell is chosen with the same as Ref. , which is a reasonable scheme to attain the helical structures in the nanowire. To characterize the thermal behavior of nanowires, we monitor the root-mean-square (rms) fluctuation of the interatomic bond distances $`\delta `$ defined by

$$\delta =\frac{2}{n(n1)}\underset{i<j}{\overset{n}{}}\frac{\left(r_{ij}^2_tr_{ij}_t^2\right)^{1/2}}{r_{ij}_t}$$
(1)
and the heat capacity $`C`$ per atom, which is related to the energy $`E`$ fluctuation, by the relation

$$C=\frac{(E^2E^2)}{nk_BT^2}$$
(2)
where $`r(i,j)`$ denotes the distance between the nuclei $`i`$ and $`j`$, $`n`$ is the total number of the atoms in the nanowire, $`k_B`$ is the Boltzman constant, $`\mathrm{}`$ indicates the thermal statistical averages in the canonical ensemble after equilibration.

The constant temperature molecular dynamics (MD) method of Nose  is employed to exploit the thermal properties of gold nanowires. The MD time step is chosen as 2.15$`\times 10^{15}`$s. The initial 105 steps are used to bring the system into equilibration and further 106 steps are used to record the thermal average of the physical quantities. We study five representative helical structures of gold nanowires with up to triple shells: single atom centered, double-chain, trigonal, tetragonal and six atoms parallelogram centered noted as S1, S2, S3, S4, S6 in Fig.1. Table I presents the starting melting temperature $`T_{ini}`$ and the overall melting temperature $`T_m`$ of these structures. The overall melting temperature $`T_m`$ reflects the stability of the whole system, while the starting melting temperature $`T_{ini}`$ describes the early stage of melting and sensitively depends on wire structures. The overall melting temperature $`T_m`$ is estimated from the curves of rms bond length fluctuation $`\delta `$, heat capacity $`C`$ and binding energy $`E`$ as functions of temperature.

As shown in Table I, the overall melting temperatures $`T_m`$ for nanowires are almost the same for all the wires and are slightly lower than the bulk melting point. Experimentally, Liu et al.  have reported the melting point of Pt nanowire is about 400oC. Lee et al. have also found the melting temperature of 4.6nm Pd nanowire is just 300oC, which is much lower than the bulk value (1445oC). Similar phenomenon is also found in the case of metal clusters. The depression of melting point can be attributed to the low dimension and large surface-to-volume ratio in these nanostructures. However, the melting temperatures obtained for gold nanowires are much higher than those of gold nanoclusters . This may be understood by the tightened helical structures in nanowires. On the other hand, the starting melting temperatures $`T_{ini}`$ are different for wires with different interior structures (see Table I). The wire S1 (as described by multi-walled structural index:18-12-6-1) starts to melt at a rather lower temperature than the other ones, while the wire S4 (21-15-9-4) has relatively higher starting melting temperature. These indicate that the different starting temperature $`T_{ini}`$ for different nanowires are related to their different interior structures. Similar overall stability for all the wires may come from their common multi-walled helical packing, despite different interior structure. Moreover, since the $`T_{ini}`$ is related to the stability of interior atomss, the rather lower starting melting temperature $`T_{ini}`$ implies that the interior melting behavior happens during melting process of the helical gold nanowires.

We discuss the structural evolution of gold nanowires during the melting process. Fig.2 gives several snap shots taken from the structural trajectories of the wire S1 at different temperatures. It is interesting to note that the interior atoms diffuse along the wire axis direction at a rather low temperature. As shown in Fig.2, the center atoms firstly move along the wire at 300K (Fig.2a). The helical structure of the outmost shell is almost invariant. With the raise of temperature, the center atoms continue moving away from the wire. The atoms in the first shell (from interior to outer) thus have fewer nearest neighbors and begin to diffuse along axis direction (Fig.2b). Similarly, the atoms in the second shell are also involved into the migration at higher temperature $`T=900`$K (Fig.2c). However, the helical structure of surface shell still exists. When the temperature is high enough, the surface atoms in the outmost shell also start to migrate and the helical structure of surface is broken at 1100K, leading to the overall melting (Fig.2d). Therefore, we conclude that interior atoms diffuse prior to surface atoms and no surface melting takes place before the overall melting in gold nanowires.

To further illustrate the argument of the interior melting and distinguish the role of surface and core atoms on the melting behavior, Fig.3 plots the rms bond length fluctuation of the surface atoms ($`\delta _s`$), core atoms ($`\delta _c`$) and the total atoms ($`\delta `$) of S1 wire as functions of temperature. Obviously, the rms bond length fluctuation of the core atoms have similar trends to that of the whole wire, but dramatic different from the surface. In the temperature range of $`3501000`$K, the $`\delta _s`$ of the surface atoms is very small and almost invariable, while $`\delta _c`$ and $`\delta `$ have substantial fluctuations. In other words, the core atoms begin to diffuse along axis and become ‘wet’, while the surface atoms remain ‘solid-like’ at this melting region. These results indicate that the melting mainly comes from the diffusion of core atoms and no surface melting happens at the beginning of the melting. Moreover, the $`\delta _c`$ fluctuates around 0.12 at low temperature, consistent with the Lindemann criterion for equilibrium melting of simple crystals, which shows that the interior melting takes place. For all the three cases, there is a rapid raise of rms bond length fluctuation in a narrow temperature region (1000-1150K), indicating that the surface atoms are involved into the melting process. Afterwards, the surface atoms play an important role on the melting of the wire at the high temperature region. Above 1150K, all of three quantities have a large, constant and smooth variation, corresponding to the completely melting status. Together with Fig.2d, we hold the surface melting means the overall melting in the multi-walled ultrathin helical nanowires. The interior melting is in preference to the surface melting and the diffusion of core atoms has a dominant effect on the melting at low temperature. The surface atoms affect the melting at high temperature and surface melting represents the overall melting.

The above interior melting behavior is obviously different from that of the bulk, clusters and crystalline nanowires, where the surface melting usually occurs before the total melting. In the bulk, after the melting of surface layers, the rest core atoms can still be seen as an analogous bulk. Much more energy is needed to make the ‘rest bulk’ molten. Therefore, the surface melting temperature can be lower than the overall melting point. In the cases of clusters and crystalline nanowires, surface atoms have fewer nearest neighbors and are thus weaker bound and less constrained compared to the core atoms. Thus, surface atoms are easy to diffuse and become liquid-like. However, the present structures of gold nanowires are helical multi-walled cylindrical. The interaction among the same shell atoms is stronger than that of neighboring shell atoms. Furthermore, to reach a well-tightened helical multi-walled structure, the helical match may cause the fewer number of atoms in the interior shells, especially for the center chain. Thus, interior atoms may have less coordinate numbers and larger interatomic distances than those of surface atoms. Therefore, core atoms in nanowires can break away from the binding sites in the wire prior to surface atoms and diffuse at a rather low temperature. In addition, in th optimization process of ground state structure, we have found that the formation of helical structure in the outer shell is earlier than the interior shell. This implies that the surface is dynamically more stable than the interior part. Experimentally, Wu et al. have found that the melting of Ge nanowire starts from the two ends of the wire and move towards the middle. These imply that the melting behavior for nanowires dramatically differs from those for nanorods and bulk.

To further clarify our ideas, we separate the function of surface atoms and core atoms by fixing the surface atoms artificially and allowing the interior atoms to move or vice verse. All the rms bond length fluctuation $`\delta `$, $`\delta _{fc}`$, $`\delta _{fs}`$ are calculated for the entire system, representing no constrained, fixing core atoms and fixing surface, respectively. We still take the wire S1 as an example. As shown in Fig.4, regardless the large difference among $`\delta `$, $`\delta _{fc}`$, $`\delta _{fs}`$ curves, the full melting temperatures (which corresponds to a large, smooth and constant rms bond length fluctuation) in the three cases are about 1150K. The $`\delta _{fs}`$ has considerable jump at low temperature (400K), while the $`\delta _{fc}`$ is almost invariable up to 800K. This indicates that surface atoms with the helical structure are thermodynamically more stable. For the case of fixing the surface and no constrained, although the $`\delta _{fs}`$ is just nearly half of the $`\delta `$, their general trends are similar. The small absolute value of $`\delta _{fs}`$ comes from the fixing surface atoms, which contributes nearly zero to the rms bond length fluctuation. Therefore, we propose there should be no essential difference between the first two cases. These results also support the melting comes from the interior atoms at low temperature and the surface melting represents the overall melting.

To check the validity of the current results, we have used the same MD code to study the melting of clusters and obtained results similar to previous Monte Carlo (MC) simulations. We have also exploited the melting process of gold nanowires by MC and observed similar interior diffuse behavior. Moreover, to examine the effect of the periodic bound condition on the melting behavior, we rescale the supercell length to 2-times and 3-times of original one. The melting temperatures are also about 1100K and the observed melting process is similar to the above results. It proves that the periodic bound condition has little effect on the simulation results. It should be further pointed out that we just limit our discussion to the helical structure in this paper. As mentioned above, this kind of helical structure is prevalent in metal nanowires in the small diameter range . Therefore, the current results on interior melting behavior in ultrathin gold nanowires is significant and might be a common feature in this kind of ultrathin metallic nanowires.

In summary, the thermal behavior of helical multi-walled gold nanowires has been studied and the main points are made as following. (1) The melting process starts from the interior region and no surface melting happens at lower temperature. We further argued that interior melting behavior happens prior to surface melting and surface melting represents the overall melting in ultrathin metallic nanowires. (2) The overall melting temperature of gold nanowires is lower than the one for the bulk, but higher than gold nanoclusters. (3) The surface and core atoms play different role in the melting behavior of these nanowires. The core atoms have a dominating effect on the melting at the beginning stage and the surface atoms are involved in the melting at higher temperature region. The core melting is closely related to the interior atomic structural characters. (4) The novel interior melting behavior is ultimately attributed to the helical multi-walled structure.

This work is financially supported by the National Natural Science Foundation of China(No.29890210, 10023001) and One-hundred-person project of Chinese Science Academic in China(2000).


text:
# 1 The model

## 1 The model

We study a two-dimensional random resistor network of total resistance $`R`$, made of $`N_{tot}`$ resistors, each of resistance $`r_n`$, disposed on a square lattice. We take a square geometry, $`N\times N`$, where $`N`$ determines the linear size of the lattice. For the comparison with resistance measurements of thin films, the value of $`N`$ can be related to the ratio between the size of the sample and the average size of the grains composing the sample. An external bias represented by a constant current $`I`$, is applied to the RRN through electrical contacts realized by perfectly conducting bars at the left and right hand sides of the network. A current $`i_n`$ is then flowing through each resistor. The RRN interacts with a thermal bath at temperature $`T_0`$ and the resistances $`r_n`$ depend linearly on the local temperatures, $`T_n`$, as:

$$r_n(T_n)=r_0[1+\alpha (T_nT_0)]$$
(1)
In this expression $`\alpha `$ is the temperature coefficient of the resistance, and $`T_n`$ is calculated by adopting a biased percolation model  as:

$$T_n=T_0+A\left[r_ni_n^2+\frac{B}{N_{neig}}\underset{l=1}{\overset{N_{neig}}{}}\left(r_li_l^2r_ni_n^2\right)\right]$$
(2)
Here, $`N_{neig}`$ is the number of first neighbours around the nth resistor, the parameter $`A`$ describes the heat coupling of each resistor with the thermal bath and it determines the importance of Joule heating effects. The parameter $`B`$ is taken to be equal to $`3/4`$ to provide a uniform heating for the case of perfect network. In the initial state (corresponding to the perfect network with no heating) all the resistors are identical: $`r_nr_0`$. Accordingly, the initial value of the resistance is $`R_0=r_0[N/(N+1)]`$. Now, we assume that two competing percolative processes act to determine the RRN evolution. The first process consists of generating fully insulating defects (resistors with very high resistance, i.e. broken resistors) with probability $`W_{D,n}=exp[E_D/K_BT_n]`$, where $`E_D`$ is an activation energy characteristic of the defect and $`K_B`$ the Boltzmann constant . The second process consists of recovering the insulating defects with probability $`W_{R,n}=exp[E_R/K_BT_n]`$, where $`E_R`$ is an activation energy characteristic of this second process. For $`A0`$, Eq. (2) implies that both the processes (defect generation and defect recovery) are correlated percolations. Indeed, the probability of breaking (recovering) a resistor is higher in the so called “hot spots” of the RRN. On the other hand, for $`A=0`$ Eq. (2) yields $`T_nT_0`$, which corresponds to random percolations . The same is true for vanishing small bias values, when Joule heating effects are negligible. As a result of the competition between these two percolations the RRN reaches a steady state or fails completely with a critical fraction of defect $`p_c`$, corresponding to the percolation threshold . In the first case, the network resistance fluctuates around an average value $`<R>`$, while in the second case, an irreversible breakdown occurs, i.e. $`R`$ diverges due to the existence of at least one continuous path of defects between the upper and lower sides of the network . We note that in the limit of a vanishing bias (random percolation) and infinite lattices ($`N\mathrm{}`$), the expression: $`E_R<E_D+K_BTln[1+exp(E_D/K_BT)]`$ provides a sufficient condition for the existence of a steady state , while it represents only a necessary condition in the case of biased percolation.

The evolution of the RRN is obtained by MC simulations carried out according to the following iterative procedure. (i) Starting from the perfect lattice with given local currents, the local temperatures $`T_n`$ are calculated according to Eq. (2); (ii) the defects are generated with probability $`W_{D,n}`$ and the resistances of the unbroken resistors are changed as specified by Eq. (1); (iii) the currents $`i_n`$ are calculated by solving Kirchhoff’s loop equations by the Gauss elimination method and the local temperatures are updated; (iv) the defects are recovered with probability $`W_{R,n}`$ and the temperature dependence of unbroken resistors is again accounted for; (v) $`R`$, $`i_n`$ and $`T_n`$ are finally calculated and the procedure is iterated from (ii) until one of the two following possibilities is achieved. In the first, electrical breakdown occurs (in numerical calculations we stop the iteration when $`R`$ increases over a factor of $`10^3`$ with respect to its initial value). In the second, the RRN reaches a steady state; in this case the iteration runs long enough to allow a fluctuation analysis to be carried out. Each iteration step can be associated with an elementary time step on an appropriate time scale (to be calibrated with experiments). In the simulations, as reasonable values of the parameters, we have taken: $`N=75`$, $`r_0=1(\mathrm{\Omega })`$, $`\alpha =10^3`$ (K<sup>-1</sup>), $`A=5\times 10^5`$ (K/W) and $`T_0=300(K)`$ if not stated otherwise. Furthermore, $`E_D=0.17(eV)`$ and values of $`E_R`$ ranging from $`0.026(eV)`$ to $`0.155(eV)`$ are used. The values of the external bias range from $`10^4I2.8`$ (A).

## 2 Results

Figure 1 reports the mean value of the RRN resistance as a function of the external current for the values of $`E_R`$ reported in the figure. Each value of $`<R>`$ has been calculated by consider the time average on a single realization and then averaging over 20 realizations. The general behavior of each set of data shows a systematic increase of the resistance up to the breakdown of the network which occurs for currents greater than the breakdown value $`I_b`$ corresponding to the last point reported in the figure. For each set of data, at increasing current values, we identify three main transport regimes, namely: a linear regime at the lowest currents, a nonlinear regime at intermediate current values, and a pre-breakdown regime at the highest currents. The linear regime corresponds to a resistance which is independent of current (Ohmic response) and extends wider the lower the recovery energy. The nonlinear regime can be described by a quadratic dependence of the resistance in an appropriate range of current values starting from a threshold value $`I_0`$ and will be analyzed with details in the following. The pre-breakdown regime covers a narrow region just before the breakdown of the network which will be better analyzed on the basis of the resistance fluctuations. The fact that by increasing the recovery energy the breakdown occurs at lower currents is explained by the concomitant growth of the average number of defects, and thus of the disorder, inside the network which implies its lower stability against breakdown.

The similar feature of the data reported in Fig. 1 suggests the possibility that all the curves belong to the same universality class. We define the value of the current $`I_0`$ by adopting the criterium of the departure of the resistance from the linear response value of $`0.05\%`$. Accordingly, Fig. 1 has been redrawn on a normalized current scale $`I/I_0`$ as is reported in Fig. 2. This figure, by showing that all curves collapse into a single one within an uncertainty of a few percent, clearly confirms the above suggestion. Furthermore, to obtain a better insight of the quadratic regime, in the insert we report the relative variation of resistance $`\left(RR_0\right)/R_0`$ as a function of the normalized current. Here, we have found that the power law:

$$\frac{R(I)}{R_0}=g(\frac{I}{I_0}),\mathrm{with}g(\frac{I}{I_0})=1+\left(\frac{I}{I_0}\right)^2,$$
(3)
is well satisfied over about two decades of current values. The evidence of such a quadratic regime is in good agreement with recent experimental results  concerning the electrical failure of carbon high-density polyethylene composites (HDPE). In particular, their failure transition is found to be of the first order which is a feature we have checked in the present model. From the study of HDPE samples with different concentration of conductive component (carbon black) it was discovered that the ratio $`R_b/R_0`$ between the breakdown resistance $`R_b`$ and the Ohmic value $`R_0`$ is a constant independent of the value of the carbon fraction and the geometry of the sample, depending only on the nature of the conducting component. This result emerges, in our context, as a consequence of the independence of the ratio $`I/I_0`$ on the initial resistance  and of Eq.(3). The electrical properties of these composites are thus nicely described by a RRN modelling.

An important source of information concerning the physical properties of the stationary regime comes from the study of the resistance fluctuations. To this purpose we analyze in the following the variance of resistance fluctuations as obtained by the simulations. Figure 3 reports the relative variance of resistance fluctuations $`\mathrm{\Delta }R^2/R^2\mathrm{\Sigma }`$, with $`\mathrm{\Delta }R^2=`$ $`R^2R^2`$, as a function of the external current for the different values of recovery energy reported in the figure. The general behavior of different sets are similar to those reported in Fig. 1 for the resistance. Also in this case we find a systematic increase of the variance at increasing current with the remarkable novelty that now the increase can be of several order of magnitude for the most ordered network. By increasing the value of $`E_R`$, the curves start with higher values, in other words, the variance grows with $`E_R`$. Each curve is initially independent of current, which corresponds to the linear regime where the variance is an intrinsic property of the RRN. Then, the variance starts increasing with current according with the establishment of the non-linear regime driven by the external current and accounted for by the bias percolation. Finally, the breakdown of the network is achieved but in a much sharp way than for the case of the resistance. We note, that the results at the two highest values of the recovery energy (open circles and full squares in the figure) are noisy already in the linear regime due to an initial network which is significantly disordered and thus less stable than those corresponding to low values of $`E_R`$.

The linear regime, occurring for small currents, when Joule heating effects are negligeable, practically corresponds to the competition of two random percolations . We have thus investigated what kind of noise characterizes this random regime. Figure 4 reports the relative variance of resistance fluctuations $`\mathrm{\Sigma }`$, as a function of the network defectiveness monitored by $`p_cp`$ where $`p_c=0.5`$ is the critical value of the fraction of defects (i.e. the percolation threshold) for a square lattice and $`p`$ the average fraction of defects under stationary conditions. The plotted data have been calculated for bias $`I<I_0`$ and each point refers to the different value of $`E_R`$ reported in the figure. For recovery energies smaller than about $`0.103(eV)`$ we obtain a quasi-perfect lattice, i.e. $`pp_c`$, while for larger values, the fraction of defects become significant and $`\mathrm{\Sigma }`$ grows following the scaling law: $`\mathrm{\Sigma }(p_cp)^k`$. Our simulations give $`k=3.2\pm 0.1`$, thus confirming previous findings obtained by taking the thermal coupling parameter $`A=0`$ . We note, that in literature  the value $`1.12\pm 0.02`$ has been found for the case of a flicker noise spectrum in disordered network subjected to random percolation. Therefore, the present noise should be interpreted as coming from the intrinsic source of fluctuations and with the property of being much more sensitive to the presence of disorder.

As for the case of the resistance, the similar features of the data reported in Fig. 3 suggest the possibility that all the curves belong to the same universality class. To this purpose, Fig. 5 reports on a log-log scale the quantity $`\mathrm{\Sigma }/\mathrm{\Sigma }_01`$ as a function of the normalized current $`I/I_0`$. Here $`\mathrm{\Sigma }_0`$ indicates the variance of resistance fluctuations for the Ohmic regime. The data are the same of those in Fig. 3 and all the curves are found to collapse onto a single one. Looking in detail the data in Fig. 5, we find that beside the quadratic regime, higher order terms are present . Indeed, by analogy with the case of resistance, the quadratic regime is identified at the lowest current ($`I<10I_0`$) as emphasized by the continuous line. At high currents, ($`I>10I_0`$) a $`4`$-th order and eventually a $`6`$-th order regime can be identified. The presence of these high order terms can be accounted for by a unique law which expands the variance of resistance fluctuations in an even power series of the external current. In particular, for large values of $`E_R`$ the breakdown happens for small $`I_b`$, so only the quadratic term of the perturbative expansion, $`\left(I/I_0\right)^2`$, is appreciable. By contrast, for small values of $`E_R`$ the breakdown happens for large $`I_b`$ and terms of order higher than the second are detected. We conclude that the noise is more sensitive than resistance to probe the pre-breakdown region, since fluctuations are here dramatically amplified with respect to negligible deviations of the resistance from the quadratic regime.

## 3 Conclusions

We have studied the steady state of a two-dimensional random resistor network resulting from the simultaneous evolutions of two competing percolations biased from an external current. The two percolations consist in a defect generation process and in a defect recovery one. The state of the network is investigated by Monte Carlo simulations over the full range of the applied stress. In particular, we have analyzed the behavior of the average resistance and of the relative variance of resistance fluctuations as a function of the external current and for different values of the recovery activation energy. We have found that both these quantities follow a scaling relation in terms of $`I/I_0`$, i.e. the ratio between the applied current and the current value corresponding to the onset of nonlinearity. The scaling exponents are found to be independent of the activation energy value, implying that all the failure processes belong to the same universality class. In the linear regime, the relative variance of resistance fluctuations has been analysed with respect to the average fraction of defects in the network and, at increasing disorder, a scaling law as a function of $`p_cp`$ has been obtained. In the nonlinear regime, the relative variance of resistance fluctuations exhibits a quadratic dependence on the current while the presence of higher order terms characterizes the pre-breakdown region. The analysis of these higher order terms is left as a subject of further investigation.

Acknowledgements

This research is performed within the STATE project of INFM. Partial support is also provided by the MADNESS II project of the Italian National Research Council and the ASI project N. 253.


text:
# Polymer collapse in the presence of hydrodynamic interactions

## I Introduction

When a polymer is placed in solution hydrophobic interactions between monomers and solvent molecules can cause it to undergo a collapse transition to a compact state. The statistical physics of the polymer transition from the extended to the collapsed state is well understood. However the dynamics of the transition and in particular the effect on that dynamics of the hydrodynamic properties of the solvent remains unclear. Therefore in this Letter we investigate numerically the dynamical behaviour of a polymer chain collapsing in a dilute solution. The collapse is measured with and without the presence of hydrodynamic interactions thus allowing a direct investigation of their effect. We find that hydrodynamics accelerates the polymer collapse. It also alters the folding pathway, allowing the folding to occur more homogeneously along the polymer chain rather than initially at the chain ends.

With the introduction of Zimm’s model it became apparent that hydrodynamics play a central role in the dynamics of polymers in dilute solution. However, understanding such interactions is difficult, analytically because they present a complicated many-body problem, and numerically because they develop on time-scales long compared to the thermal fluctuations of the monomers.

Theoretical work on the dynamics of polymer collapse can be divided into two approaches. Phenomenological models balance the driving and dissipative forces to give scaling laws. They involve assumptions about how the collapsed state develops on which there is no consensus. Several authors have considered models which are based on a solution of the Langevin equation. Of particular interest is work by Pitard and by Kuznetsov et al who find the inclusion of hydrodynamics, modelled by a preaveraged Oseen tensor, speeds up the collapse.

Simulations on polymer collapse, using Monte Carlo or Langevin approaches, have not included the hydrodynamic effects of the solvent. Very recent work has shown that it is now possible to use molecular dynamics simulations with an explicit solvent to model the collapse transition if powerful computational resources are available. In an interesting recent paper Chang and Yethiraj compared molecular dynamics simulations of a polymer in a solvent to Brownian dynamics simulations. They attempted to match the parameters in the two simulations and hence to compare collapse with and without hydrodynamics. Here we use a hybrid approach where the solvent is modelled by a Malevanets-Kapral method and the polymer by molecular dynamics in our investigation of the hydrodynamics of polymer collapse.

## II Simulation details

Modelling a dilute polymer solution is a difficult task because of the existence of widely differing time scales. The dynamical properties of polymers can be dominated by hydrodynamic interactions between different parts of the polymer chains. In contrast with the time scale of thermal fluctuations of individual monomers, these interactions are long-ranged and evolve slowly. Therefore it is computationally too expensive to reach hydrodynamic time scales using molecular dynamics simulations for both the polymer and the solvent molecules.

To overcome this problem we use a hybrid simulation approach where the equations of motion of the polymer alone are solved using a molecular dynamics algorithm. The solvent is modelled using a mesoscale approach, developed by Malevanets and Kapral. This ignores the molecular detail of the solvent but preserves its ability to transmit hydrodynamic forces. The polymer can be thought of as moving within a “hydrodynamic heat bath”.

The polymer chain is modelled by beads connected via non-harmonic springs with adjacent beads along the chain backbone representing an effective Kuhn length of the polymer chain. These finitely extensible springs are represented by the FENE potential

$`V_{FENE}\left(r\right)={\displaystyle \frac{\kappa }{2}}R_{0}^{}{}_{}{}^{2}\mathrm{ln}\left[1\left({\displaystyle \frac{r}{R_0}}\right)^2\right],r<R_o.`$ (1)
A Lennard-Jones potential which acts between all the polymer beads is used to model the excluded volume of the monomers and a long range attraction which drives polymer collapse

$`V_{LJ}\left(r\right)=4\epsilon \left[\left({\displaystyle \frac{\sigma }{r}}\right)^{12}\left({\displaystyle \frac{\sigma }{r}}\right)^6\right].`$ (2)
We take $`\epsilon =1.0`$, $`\sigma =1.0,\kappa =30`$ and $`R_0=2`$ where parameters and results are quoted in reduced Lennard-Jones units.

Newton’s equations of motion for the polymer are integrated using the time reversible velocity Verlet algorithm. The molecular dynamics time step is chosen to be $`\delta t=0.002t_s`$ where $`t_s`$ is the interval between solvent collision steps, defined below.

The solvent is modelled by a large number $`N=131072`$ of point-like particles which move in continuous space with continuous velocities but discretely in time . The algorithm is separated into two stages. In the first of these, a free streaming step, the positions of the solvent particles at time $`t`$, $`𝒙_i(t)`$, are updated simultaneously according to

$`𝒙_i(t+t_s)=𝒙_i(t)+𝒗_i(t)t_s`$ (3)
where $`𝒗_i(t)`$ is the velocity of a particle.

The second component of the algorithm is a collision step which is executed on both solvent particles and polymer beads. The system is coarse-grained into $`L^3`$ unit cells of a regular cubic lattice. In this simulation $`L=32`$ is used. There is no restriction on the total number of solvent or polymer particles in each cell, although the total number of particles is conserved. Multiparticle collisions are performed within each individual cell of the coarse-grained system by rotating the velocity of each particle relative to the centre of mass velocity $`𝒗_{cm}(t)`$ of all the particles within that cell

$`𝒗_i(t+t_s)=𝒗_{cm}(t)+𝑹\left(𝒗_i(t)𝒗_{cm}(t)\right).`$ (4)
$`𝑹`$ is a rotation matrix which rotates velocities by $`\theta `$ around an axis generated randomly for each cell and at each time step. In the present calculations we take $`\theta =\frac{\pi }{2}`$.

Note that the collision step preserves the position of the solvent and polymer beads. It transfers momentum between the particles within a given cell while conserving the total momentum of these particles. Because both momentum and energy are conserved locally the thermohydrodynamic equations of motion are captured in the continuum limit. Hence hydrodynamic interactions can be propagated by the solvent and, because the polymer beads are involved in the collisions, to the polymer. Note, however, that molecular details of the solvent are excluded: this allows the hydrodynamic interactions to be modelled with minimal computational expense.

The volume in phase space is invariant under both the free streaming and collision steps. Hence the system is described by a microcanonical distribution at equilibrium. The initial solvent distribution was generated by assigning positions randomly within the system with an average density $`\rho =4`$ particles per unit cell. The velocities were assigned from a uniform distribution which relaxed rapidly ($`t_s100`$) to the equilibrium Maxwell-Boltzmann form.

A particularly useful feature of the Malevanets-Kapral algorithm is the ease with which hydrodynamic interactions can be “turned-off” thus replacing the hydrodynamic heat bath by a Brownian (random) heat bath. This is achieved by randomly interchanging the velocities of all the solvent particles after each collision step, thus relaxing the constraint of local momentum conservation to a global one. Accordingly the velocity correlations which result in hydrodynamic interactions disappear from the fluid. Running simulations with the same initial conditions and parameter values, but with hydrodynamics present or absent, greatly faciliates pinpointing the effect of the hydrodynamic interactions.

## III Results

Our aim is to study the dynamics of polymer collapse from an extended to a compact state. The collapse trasnsition is driven by the attractive Lennard-Jones interactions between the polymer beads. The transition is rounded and shifted by the finite chain length, and takes place at $`k_BT1.8`$ for chains of length $`N_p=100`$.

Initial polymer configurations were chosen at random from long runs on a polymer chain at equilibrium in a solvent with $`k_BT=4`$. Each of these extended configurations were placed in a solvent at equilibrium at $`k_BT=0.8`$. This value ensures collapse but prevents the chain collapsing so quickly that hydrodynamic interactions do not have sufficient time to develop. The rate of polymer collapse was measured by monitoring the variation of the radius of gyration $`R_g`$ of the chain with time

$`R_{g}^{}{}_{}{}^{2}(t)`$ $`=`$ $`{\displaystyle \frac{1}{N_p}}{\displaystyle \underset{i=1}{\overset{N_p}{}}}\left(R_i(t)R_{cm}(t)\right)^2,`$ (5)
$`R_{cm}(t)`$ $`=`$ $`{\displaystyle \frac{1}{N_p}}{\displaystyle \underset{i=1}{\overset{N_p}{}}}R_i(t).`$ (6)
A typical numerical result is shown in Figure 1 for a chain of length $`N_p=100`$ with and without hydrodynamics.

The collapse time was estimated as the time when the equilibrium radius of gyration was first attained. Results for $`\tau `$ are shown for chains of varying lengths in Table 1. $`\tau `$ was averaged over 20 initial configurations for $`N_p=40`$, $`60`$ and $`100`$, and 5 for $`N_p=200`$, and the collapse time for hydrodynamic ($`\tau _H`$) and Brownian ($`\tau _B`$) heat baths compared. Variations in collapse times between different runs are large as expected. However, it is strikingly apparent that hydrodynamic interactions speed up the rate of collapse by a factor$`2`$ for each polymer length. For all but 6 of the 65 runs and for all the $`N_p=100`$ and $`200`$ runs the collapse was faster with hydrodynamics switched on. The results in Table 1 also illustrate that the collapse time increases with chain length as expected ($`\tau _BN_p^{1.6\pm 0.2}`$, $`\tau _HN_p^{1.2\pm 0.3}`$). Values predicted in the literature for the exponents relating the collapse time to chain length vary widely. It is far from obvious that these exponents are universal: they appear to depend on quench depth and the details of the model. Moreover the final radius of gyration was recorded for both the Brownian and hydrodynamic collapses and found to agree as the equilibrium polymer properties should be independent of the nature of the heat bath.

Figure 2 compares typical collapse pathways with and without hydrodynamics for $`N_p=100`$ and a final temperature $`k_BT=0.8`$. In this and the majority of other runs a qualitative difference was observed in the collapse mechanism. In the Brownian solvent the ends of the polymer tend to collapse first forming a dumbbell shape. With hydrodynamics the collapse takes place more evenly along the chain. We speculate that this occurs due to ’slip-streaming’. The movement of beads towards each other, which is initiated by the attractive Lennard-Jones potential, is enhanced by hydrodynamic interactions. Once a bead starts moving in a given direction and locally drags fluid with it, it is easier for neighbouring beads to move in the same direction. Similar behaviour albeit on much larger length scales has been observed experimentally by Bartlett et al for a pair of colloidal particles in solution.

## IV Discussion

We have shown that it is possible to measure the dynamics of the collapse transition of a polymer in a solvent using a hybrid mesoscale/molecular dynamics algorithm. Hence it has been possible to show that hydrodynamic interactions speed up the collapse of the polymer chain. The hydrodynamics alters the collapse pathway, allowing the folding to occur more homogeneously along the chain, rather than initially at the chain ends.

Qualitative features of the collapse pathways are in agreement with those reported by Chang and Yethiraj who recently compared molecular dynamics simulations of a polymer in a solvent, which include hydrodynamics, to Brownian dynamics simulations, which did not. These authors also found that the Brownian simulations could become trapped in a metastable free energy minimum. It would be interesting to address whether this is due to deeper quenches than those considered here or to the different simulation approaches employed.

The solvent is modelled using the Malevanets-Kapral algorithm which sustains hydrodynamic modes but does not include molecular interactions. Hence we caution that we cannot investigate the late stages of collapse where trapping of the water molecules in the polymer chain may be important. Moreover collapse is driven by attractive Lennard-Jones interactions between the polymer beads rather than repulsive monomer-solvent interactions. It would be of interest to include the solvent near the polymer in the more realistic molecular dynamics updating and the possibility of doing this is under investigation.

Finally it may be of interest to consider whether the early stages of the folding of proteins in aqueous solution might be influenced by the hydrodynamic interactions of the solvent increasing the cooperativity of the collapse.

Acknowledgements We thank A. Malevanets, C. Pooley and P. Warren for helpful discussions.