text:
# Single crystal growth and investigation of 𝐍𝐚_𝐱𝐂𝐨𝐎_𝟐 and 𝐍𝐚_𝐱𝐂𝐨𝐎_𝟐𝐲𝐇_𝟐𝐎

## I Introduction

The recently discovered superconductor $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2\mathrm{yH}_2\mathrm{O}`$ has attracted considerable attention as being the first superconducting cobaltite and due to evidence for important electronic correlations 1 . Its superconducting transition temperature of maximum $`\mathrm{T}_\mathrm{c}`$5 K exhibits a composition dependence, decreasing for both underdoped and overdoped materials, as observed in the cuprates. This suggests that a detailed characterization of the electronic and magnetic behavior of this new type of materials and their interplay with structural peculiarities may contribute to a more fundamental understanding of the high $`\mathrm{T}_\mathrm{c}`$ superconductivity in cuprates, such as the layered crystal structure 2 , magnetic susceptibility 3 ; 4 , electronic anisotropy 5 , irreversibility lines and superconducting parameters $`\lambda _{\mathrm{ab}}`$, $`\mathrm{H}_{\mathrm{c2}}`$.

The physical properties are determined by the amount of sodium and water in $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2\mathrm{yH}_2\mathrm{O}`$ via their influence on the Co valence state and the resulting local distortions of the $`\mathrm{CoO}_6`$ octahedra present in the structure. The compound $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ shows manifestations of frustration in its physical properties, because Co occupies a triangular lattice and the exchange interactions are antiferromagnetic 6 ; 7 . The close relationship between structural and electronic properties establishes itself also on the Fermi-Surface that is expected to show features of strong nesting tanaka03 ; johannes04 . Furthermore, there exists evidence that superconductivity might have an unconventional order parameter tanaka03 . When water is introduced into the compound, it exhibits chemical and structural instabilities. The $`\mathrm{Na}^+`$ ions are at least partially mobile at elevated temperatures and may order at low temperatures in dependence of composition and sample treatment. The same is proposed for the $`\mathrm{Co}^{4+}/\mathrm{Co}^{3+}`$ charge state on the triangular lattice huang04 .

All of these features make the study of the physical and chemical properties of the compounds difficult. However, nearly all the studies reported so far have been made on either sintered specimens or poorly characterized polycrystalline samples. These materials are often found to be inhomogeneous, especially with respect to the distribution of $`\mathrm{Na}^+`$ ions in the lattice and in the intergranular spaces, the $`\mathrm{Na}_2\mathrm{O}`$ decomposition of the compound, and with respect to the intercalation of water. In polycrystalline samples water may in part accumulate in intergranular spaces. The possible presence of second phases leaves a degree of uncertainty when interpreting the structure and the electrical and superconducting properties as function of composition, i.e. establishing a superconducting phase diagram. This problem may be circumvented by the use of high quality single crystals whose chemical composition and crystal structure can be properly determined.

When growing $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ single crystals, considerable difficulties appear on account of the high $`\mathrm{Na}_2\mathrm{O}`$ vapor pressure, which increases exponentially from $`10^5`$ to $`10^3`$ torr with heating to temperatures between 500 and 800 C, followed by a noticeable evaporation. Therefore ceramic powders were synthesized using additional $`\mathrm{Na}_2\mathrm{CoO}_3`$ for compensation of the Na loss during heating 2 ; 8  or a ”rapid heat-up” technique 9  to avoid the formation of a non-stoichiometric compound. Solution growth by using NaCl flux 10  was performed that unfortunately lead only to thin crystals ($`<`$0.03 mm) or non-stoichiometric and possibly contaminated samples. To obtain the superconducting phase the crystal must be hydrated by chemical oxidation 1 . The compound can be partially decomposed during the oxidation process, leading to highly defective crystals containing Na-poor phases.

In this contribution we present a single crystal growth method to prepare large and high quality single crystals and demonstrate their chemical, thermal, and structural behavior and electric properties under dry and humid conditions.

## II Experimental

Single crystals were grown in an optical floating zone furnace (Crystal System Incorporation, Japan) with 4 x 300 W halogen lamps installed as infrared radiation sources. Starting feed and seed materials were prepared from $`\mathrm{Na}_2\mathrm{CO}_3`$ and $`\mathrm{Co}_3\mathrm{O}_4`$ of 99.9% purity with the nominal composition of $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$, where *x*=0.50, 0.60, 0.70, 0.75, 0.80, 0.85 and 0.90, respectively. Well-mixed powders were loaded into alumina crucibles and heated at 750 C for a day. The heated powders were reground and calcined at 850 C for another day. They were then shaped into cylindrical bars of $``$ 6x100 mm by pressing at an isostatic pressure of $``$70 Mpa and then sintered at 850 C for one day in flowing oxygen to form feed rods.

The sintered feed rod was premelted with a mirror scanning velocity of 27 mm/h by travelling the upper and lower shafts, respectively, to densify the feed rod. After premelting the $``$20 mm rod long was cut and used as a first seed and hereafter the grown crystal was used as a seed. The feed rod and the growing crystal were rotated at 15 rpm in opposite directions. In an attempt to reduce the volatilization of Na and obtain stoichiometric and large crystals, we applied travelling rates of 2 mm/h under pure oxygen flow of 200 ml/min throughout the growing procedure.

To obtain superconductivity the Na was partly extracted by placing crystals in a 6.6 mol $`\mathrm{Br}_2/\mathrm{CH}_3\mathrm{CN}`$ solution during 100 hours and then washed out the solution by $`\mathrm{H}_2\mathrm{O}`$ or $`\mathrm{D}_2\mathrm{O}`$ 1 . Alternatively, the electrochemical technique was also applied, using an aqueous solution of NaOH as an electrolyte to partially extract Na 11 . This technique needs a longer time for deintercalation and the resulting superconducting transitions are generally sharper.

Thermogravimetric and differential thermal analysis (TG-DTA) was performed to study the melting behavior as well as the time and temperature dependence of the water loss in the compound. Single crystal XRD measurements were carried out with the X-ray diffractometer (Philips PW 1710) using Cu $`K_\alpha `$ radiation, a scanning rate of 0.02 degrees per minute, and $`\theta `$ \- 2 $`\theta `$ scans from 5 to 90 degrees. The lattice parameters obtained from XRD data were refined using the commercial program PowderCell. The chemical composition was determined by energy dispersive X-ray (EDX) including microanalysis. The as-grown crystals were cleaved parallel to the *a* axis and the composition of Na and Co was determined across the cleaved section.

## III Results and discussions

### III.1 DTA-TG analysis

The melting behavior of a $`\mathrm{Na}_{0.7}\mathrm{CoO}_2`$ crystal was investigated by DTA-TG measurements and with a high temperature optical microscope. A small single crystal was loaded in a Pt crucible and then heated in the DTA-TG apparatus (NETZSCH STA-449C) at 7.5 C/min up to 1200-1250 C in flowing oxygen. Two peaks with onset $`\mathrm{T}_1`$=1035 C and $`\mathrm{T}_2`$=1092 C on the DTA curve correlate very well with a weight loss observed by TG, as shown in Figs. 1(a) and (b), respectively. No weight changes were observed below 980 C.

The high temperature optical microscope study of the crystal reveals that the liquid phase appears at $`\mathrm{T}_1`$=1035 C but the solid phase (crystal) still remains up to $`\mathrm{T}_2`$=1092 C, exhibiting a coexistence of solid and liquid phases. Our investigations indicate that thermal decomposition of $`\mathrm{Na}_{0.7}\mathrm{CoO}_2`$ is accompanied by a weight loss down to 96.2 wt.% and the compound decomposes into a sodium-rich liquid and a cobalt-rich solid phase 12 ; 13 , assuming chemical reaction of the melt as follows, $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ $`>`$ CoO + Liquid (Na-rich) + *x*O2 \- (T$`>`$1100 C) $`>`$ \[$`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_{1.65}`$\]\* (\*=homogeneous melt).

The dissolution of CoO occurs in the Na-rich melt and by taking up oxygen from the environment. Thus, it results in an increase of weight because the oxidation state of cobalt is significantly higher in the melt, i.e., $`\mathrm{Co}^{+2.7}`$, according to the chemical reaction equation shown above. A sharp weight loss down to 95.5 wt.% occurs at $`\mathrm{T}_2`$=1092 C and the melt starts to get homogeneous and stable. A nearly constant weight of 95.5 - 94.9 wt.% is obtained in the temperature range of 1120 - 1200 C, where the melt becomes homogeneous. A monotonic weight loss of melt may happen under constant heating. Evidently, the melting behavior described so far indicates that the compound melts incongruently.

### III.2 Crystal growth, morphology and composition

When sintering $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ at high temperatures the decomposition of $`\mathrm{Na}_2\mathrm{O}`$ can cause a weight loss and lead to an inhomogeneous compound. During growth a white powder of $`\mathrm{Na}_2\mathrm{O}`$ is observed to volatilize and deposit on the inner wall of the quartz tube. By weighing a total weight loss $``$6 wt.% is estimated. This value is in agreement with the data shown in Fig. 1(b), neglecting the tiny loss of oxygen caused by the change of cobalt valence state. The main loss of $`\mathrm{Na}_2\mathrm{O}`$ is found during the premelting procedure, estimated to be $``$5 wt.%. This is probably caused by incompletely reacted $`\mathrm{Na}_2\mathrm{O}`$ when the mixtures are calcined and sintered prior to premelting. An additional $``$1 wt.% loss is estimated after the crystal growth.

During growth, we observed that the molten zone is rather stable and $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ single crystals with a sodium content of *x*=0.70 - 0.90 is easily formed. In contrast, it is hard to grow the compound with *x*=0.50 - 0.60. The analysis of EDX indicates that as-grown $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ with *x*$`<`$0.6 consists of multi phases like, $`\mathrm{Na}_2\mathrm{O}`$, $`\mathrm{CoO}_2`$ and Na-poor phases.

Large and high quality single crystals were obtained with a growth rate of $``$2 mm/h in flowing oxygen atmosphere. Figure 2 shows a typical crystal ingot of $`\mathrm{Na}_{0.7}\mathrm{CoO}_2`$. The *001* face is readily cleaved mechanically due to the layered structure of the compound. Figure 3(a) displays one halve of the crystal platelets with the *001* face of several $`\mathrm{cm}^2`$ area cleaved from an ingot with a sharp scalpel. Figure 3(b) shows the other half after water intercalation. Crystal grains are found to grow preferably along the *a* crystallographic axis, parallel to the rod axis.

A convex growth interface is observed to be the boundary between the regions of smaller and larger diameter of the crystal ingot, as shown in Fig. 3. The smaller one is formed when the molten zone was at lower temperatures and the larger one is formed at higher temperatures. Therefore, the unequal diameter for an ingot indicates that temperature fluctuations occurred at the molten zone during growth. According to the temperature-composition relationship a fluctuation of heating temperature may result in a composition change of the molten zone, leading to an inhomogeneous compound. Figure 3 (a) shows that many tiny crystal grains of $`\mathrm{CoO}_2`$ gather at the boundary of the growth front. These grains can be removed after the deintercalation treatment, as shown in Figs. 3 (b) and (c), respectively. It is important to maintain a stable molten zone by applying a constant heating temperature for growing a high quality single crystal.

Single crystal XRD powder diffraction of $`\mathrm{Na}_{0.7}\mathrm{CoO}_2`$ gives the following lattice parameters and cell volume of the space group P63/mmc: *a*=2.8278(4) Å, *c*= 10.9073(0) Å, *v*=75.54 Å3. Sintered powders show a slightly different XRD patterns with *a*=2.8230(2) Å, *c*= 10.9628(8) Å 1 .

For the determination of Na/Co compositions using EDX the crystal was scanned through a segment of 3 mm along the growth direction taking the average of four measured points in the central and edge region of the crystal. It is found that the Na content varied with the temperature fluctuations during growth. At the beginning of the growth the temperature fluctuations are high and cause a high variation of the composition, $`\mathrm{\Delta }`$x$``$0.11, determined within 2 cm from the seeding part of the ingot. Following the seeding part the variation of Na content is smaller, $`\mathrm{\Delta }`$x$``$0.06, when the molten zone is maintained in a stable state by a constant temperature. A volatilized white $`\mathrm{Na}_2\mathrm{O}`$ powder accumulates on the inner wall of the quartz tube after the growth is completed. The loss of Na may result to a reduction of its concentration in the as grown crystal. The white $`\mathrm{Na}_2\mathrm{O}`$ powders also form on the surface of crystals if the samples are stored under ambient conditions. Therefore the grown crystals must be stored in an evacuated container or a desiccator to avoid decomposition.

### III.3 Sodium extraction and water intercalation

The superconducting phase $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2\mathrm{yH}_2\mathrm{O}`$ (0.26$``$ *x*$``$ 0.42, y$``$1.3) is obtained by chemically extracting (deintercalation) sodium from $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$, followed by hydration. Single crystals of $`\mathrm{Na}_{0.7}\mathrm{CoO}_2`$ were placed in the oxidizing agent $`\mathrm{Br}_2/\mathrm{CH}_3\mathrm{CN}`$ for around 100 hours, and then washed with acetonitrile. The change of sodium content of the resulting crystals is generally proportional to the bromine concentration in the $`\mathrm{CH}_3\mathrm{CN}`$ agent 3 . A crystal of $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ with *x*=0.3 can result from a 6.6 mol $`\mathrm{Br}_2/\mathrm{CH}_3\mathrm{CN}`$ treatment only after a long extraction time of more than a week. The Na extraction of the compound was also carried out by the electrochemical method. The composition $`\mathrm{Na}_{0.3}\mathrm{CoO}_2`$ could be achieved in aqueous solution of NaOH using a constant current of 0.5 mA and a voltage of 1.0 V for over 10 days. Comparing the published 1 to 5 days extraction time for the $`\mathrm{Na}_{0.7}\mathrm{CoO}_2`$ powders 1 ; 3  the crystal needs longer time to complete the extraction, also depending on the dimension of the sample. The slower extraction reaction in single crystals compared to a powder is related to the higher perfection and longer period of Na-O-Na within the NaO layers. In powder samples variations of the bonding are expected on a nanometer scale. Before and after the extraction treatment the sodium composition distribution across the crystals was determined by EDX.

The superconducting phase is then obtained by immersing the deintercalated samples for a week in $`\mathrm{H}_2\mathrm{O}/\mathrm{D}_2\mathrm{O}`$, or sealing the samples in a water vapor saturated container at room temperature. After hydration/deuteration a large increase in thickness is visible to the naked eyes and the morphology exhibits ”booklet”-like layered cracks perpendicular to the *c* axis. The sodium intercalant layer expands with decreasing Na content, because the Na removal results in Co oxidation ($`\mathrm{Co}^{3+}`$ ions oxidized to smaller $`\mathrm{Co}^{4+}`$ ions) and thus the $`\mathrm{CoO}_2`$ layers are expected to shrink 1 ; 2 . This suggests that a decreased bonding interaction between layers with decreasing Na content may result in a readily cleaving plane. Figures 4(a) and (b) illustrate the layered structures of the non- and fully hydrated phases of the compound, respectively. According to the layered structure of $`\mathrm{Na}_\mathrm{x}\mathrm{Co}_2`$ 14 , 24 Co-O bonds in the layer of $`\mathrm{CoO}_2`$ form an octahedral $`\mathrm{CoO}_6`$ with high bonding energy and therefore the layer is structurally robust, while the bonding energy is much weaker in the NaO layer because the $`\mathrm{Na}^+`$ mobility is high and the number of bonds Na1-O is 6 and Na2-O only 4/3. Therefore the Na layer readily collapses to leave the $`\mathrm{CoO}_2`$ layer terminated as the outer-most surface of the crystal. Moreover, the insertion of water induces two additional intercalated layers of $`\mathrm{H}_2\mathrm{O}`$ between NaO and $`\mathrm{CoO}_2`$, i.e., from Co-Na-Co to $`\mathrm{Co}\mathrm{H}_2\mathrm{O}\mathrm{Na}\mathrm{H}_2\mathrm{O}\mathrm{Na}\mathrm{Co}`$, which results in the expansion of the c lattice constant from 11.2 to 19.6 Å. The hydrogen in the new intercalated $`\mathrm{H}_2\mathrm{O}`$ layer bonds extremely weakly to the NaO and CoO layers. Therefore this compound is exceptionally unstable with respect to a change of thermal, mechanical or humidity environment. Under ambient, oxidizing and humid conditions the crystal exhibits three typical morphologies, as shown in Figs. 5(a), (b) and (c), respectively.

### III.4 Behavior of the crystals upon hydration and dehydration

To study the phase formation and problems with rehydration we have performed detailed x-ray diffraction experiments as function of time. Single crystals of non-hydrated $`\mathrm{Na}_{0.3}\mathrm{CoO}_2`$ have been partially hydrated in humid air to form $`\mathrm{Na}_{0.3}\mathrm{CoO}_20.6\mathrm{H}_2\mathrm{O}`$ under ambient conditions. These crystals were immersed in water with the *c* surface exposed to air. Each measurement of the *002* reflections was carried out for a short time of $``$3 minutes to avoid that the x-ray emits off the sample surface due to the rapid expansion of crystal thickness during the hydrating process. We show in Fig. 6 that the increase of intensity of the *002* reflection (corresponding to *y*=1.3) is accompanied by a decrease of intensity of the other two *002* reflections (corresponding to *y*=0.6 and 0, respectively). This is indicative of the increase of the volume of hydrated phase *y*=1.3 while decreasing the non-hydrated and partial hydrate phases, *y*=0 and 0.6. After a one day hydration the sample becomes completely ”wet”. This is indicated by the *002* reflections for the *y*=1.3 and 0.6 phases and the disappearance of the *002* corresponding to the non-hydrated phase *y*=0. Note that the *002* reflection of the partial hydrate with *y*=0.6 still exists until a further hydration up to 10 days. In the meantime the XRD pattern shows a continuously increasing volume fraction of *y*=1.3 phase reaching approximately 100%. This suggests that the initial hydration process takes place with two-water molecules per unit cell (corresponding to *y*=0.6) inserted into the Na plane, and followed by a group of four (corresponding to *y*=1.3) located in vacancies of the extracted Na sites.

The hydrated volume of a crystal is a function of hydration time since the water diffuses gradually in the *ab* plane whereas the *c* face is robust and no diffusion path exists along the *c* axis. A fully hydrated phase *y*=1.3, exhibiting superconductivity, is achieved after 10 days for a crystal of 2x2 $`\mathrm{mm}^2`$ in the *a*/*b* direction, concomitant with the disappearance of the partial hydrates of the other two *002* reflections, as shown in Fig. 6. This indicates that Na-vacant sites can be completely occupied by water molecules when a longer diffusion time is allowed for creating a water-saturated environment during the hydration process.

A thermogravimetric analysis can be used to monitor the loss of water and stability of the hydrated phases as function of temperature and time. Figure 7 shows the temperature dependence of water loss for an over hydrated sample (*y*=1.8) under a slow heating rate of 0.3 C /min in flowing oxygen. The fastest loss of water occurs between 20 and 50 C and is then again increasing upon further heating. The sample weight is characterized by several relatively broad plateaus in weight that are attributed to phases with *y*=1.8, 0.9, 0.6, 0.3 and 0.1 for temperatures near 28, 50, 100, 200 and 300 C, respectively. These data suggest the existence of the following majority phases at the respective temperatures: $`\mathrm{Na}_{0.3}\mathrm{CoO}_21.8\mathrm{H}_2\mathrm{O}`$ (*c*=19.75 Å), $`\mathrm{Na}_{0.3}\mathrm{CoO}_20.9\mathrm{H}_2\mathrm{O}`$ (*c*=13.85 Å), $`\mathrm{Na}_{0.3}\mathrm{CoO}_20.6\mathrm{H}_2\mathrm{O}`$ (*c*=13.82 Å) and $`\mathrm{Na}_{0.3}\mathrm{CoO}_20.3\mathrm{H}_2\mathrm{O}`$ (*c*=12.49Å) and $`\mathrm{Na}_{0.3}\mathrm{CoO}_20.1\mathrm{H}_2\mathrm{O}`$ (*c*=11.20Å). A change of $`\mathrm{\Delta }`$y=0.33 corresponds to the removal of one water molecule per unit cell. A complete removal of water (*y*=0) from the compound takes place at temperatures over 600 C.

Analysis of DTA-TG confirms that the fully hydrated superconducting phase (*y*=1.3) is rather unstable between 20 and 50 C due to the rapid loss of water. However, the superconducting phase can be reformed easily by rehydration. The inset of Fig. 7 shows the time dependence of the loss of water at room temperature. A rapid loss of water from *y*=1.8 to 0.6 occurs in one and a half hours and then no further weight loss is observed during one day. XRD reveals that the flat plateau corresponding to a semi hydrated $`\mathrm{Na}_{0.3}\mathrm{CoO}_20.6\mathrm{H}_2\mathrm{O}`$ with *y*=0.6 appears to be a metastable phase under ambient conditions. This observation is in conflict with that for polycrystalline powders 2 .

The fact of the matter is that a single crystal is only one domain containing crystal-water while powder samples consist of countless tiny grains with additional intergrain water. Furthermore, the loss or absorption of water in a large single crystal is slower and more well-defined compared to powders due to the diffusion path along the Na-layer.

### III.5 Mixture of cells

As shown in the TGA data, a fully hydrated sample can release water under ambient conditions to form partial hydrates. On the other hand, moisture from air can also be absorbed by a non-hydrated sample to form hydrated phases. In order to quantitatively identify the hydrated phases in the crystal, the fully hydrated or non-hydrated samples were placed in a humid or dry atmosphere. They were then investigated by XRD and the results show a clear shift in the positions of the *00l* reflections for both samples, yielding several variations in the *c* axes of the unit cells, as shown in Figs. 8(a-c). The fully hydrated phase of $`\mathrm{Na}_{0.3}\mathrm{CoO}_21.3\mathrm{H}_2\mathrm{O}`$ is exposed to dry air at ambient conditions for 2 and 5 days, resulting in the formation of four and three mixed phases, as shown in Figs. 8(a) and (b), respectively. Exposing a non-hydrated $`\mathrm{Na}_{0.3}\mathrm{CoO}_2`$ to humid air for 5 days results in two additional phases of $`\mathrm{Na}_{0.3}\mathrm{CoO}_20.3\mathrm{H}_2\mathrm{O}`$ and $`\mathrm{Na}_{0.3}\mathrm{CoO}_20.6\mathrm{H}_2\mathrm{O}`$, as shown in Fig. 8(c). Note that the phase with *y*=0.6 is always observed in the crystals under ambient conditions. Nevertheless, to prevent a fully hydrated phase from losing water we suggest that the samples be stored in a saturated humid atmosphere. Non-hydrated samples can be stored in an evacuated container to prevent decomposition from absorbing water in air.

### III.6 Electrical properties

We have performed transport experiments of the out-of-plane resistivity $`\rho _\mathrm{c}`$ as function of time to further characterize the hydration dynamics. The initially non-hydrated single crystal $`\mathrm{Na}_{0.3}\mathrm{CoO}_2`$ was immersed in D.I. water throughout the whole experiment. In contrast to measurements of powder samples in humid air, this configuration reduces the influence of internal and external surfaces and allows a full and continuous hydration. Furthermore, the out-of-plane resistivity $`\rho _\mathrm{c}`$ measured in four-point geometry should be less affected by the intrinsic Na ion conductivity.

In Fig. 9 and its inset plateaus and steps of the resistivity are given with increasing hydration time followed by a continuous and then more gradual decrease. Noteworthy is the similar resistivity reached after 8 hours compared to the initial state. This implies that, although irreversible by nature, the hydration process has also reversible stages. During hydration the water molecules enter the compound to intercalate between the $`\mathrm{CoO}_2`$ and Na layers. The $`\mathrm{CoO}_2`$ layers are separated by a trilayer of $`\mathrm{H}_2\mathrm{O}`$/Na/$`\mathrm{H}_2\mathrm{O}`$. The dominant effect is to expand the c-axis lattice parameters, from about 11.2 Å without water (non-superconducting phase $`\mathrm{Na}_{0.3}\mathrm{CoO}_2\mathrm{yH}_2\mathrm{O}`$, *y*=0) to 12.5 Å (semi hydrated phase, *y*=0.6) and then to 19.5 Å when the system becomes superconducting (*y*=1.3) (see Fig. 6). This is related to an enlarged interlayer distance from 5.6 Å to 6.9 Å and then to 9.8 Å. The stepwise increase from $`\mathrm{t}_0`$ to $`\mathrm{t}_5`$ occurs at a constant chemical doping level and suggests that the increase of spacing between the $`\mathrm{CoO}_2`$ layers leads to the resistivity increase. We propose that the rapid steps followed by plateaus map the dynamics of domains of hydration that are pinned at the random potential of the local defect landscape. The more gradually decrease of resistivity for t$`>`$12 hours must therefore be related to a more homogeneous state of the system.

There is an interesting correlation between the $`\mathrm{CoO}_2`$ layer and the superconducting transition temperature, i.e., the Co-O distance decreases with increasing $`\mathrm{T}_\mathrm{c}`$ 14 . The layer thickness of $`\mathrm{CoO}_2`$ changes from 1.93 Å for non-superconducting phase to 1.84 Å for superconducting phase 15 . It is clear that the layer thickness of $`\mathrm{CoO}_2`$ shrinks after reaching full hydration or doping to *y*=1.3 and a change in the oxidation state of Co occurs, because the oxidation state directly corresponds to the carrier density in the $`\mathrm{CoO}_2`$ layers and some electron transfer off the cobalt 1 . Recently, also a relation between the *c* axis parameter and the superconducting transition temperature has been discussed milne04 . Nevertheless, the change of the room temperature resistivity during the hydration process is not only influenced by the expansion of spacing in the system but also by a redistribution of the charge in $`\mathrm{CoO}_2`$ layer.

### III.7 Superconductivity

The superconducting properties of the fully hydrated single crystals were characterized by a.c. susceptibility using a SQUID magnetometer. Zero field cooling measurements with 10 Oe applied field are presented in Fig. 10. Superconducting transition temperatures $`\mathrm{T}_\mathrm{c}`$=2.8 - 4.9 K with *x*$``$0.28 - 0.42 are observed, respectively. Some of our samples showed strong diamagnetic a.c. signals and a superconducting volume fraction as high as 20%. The transition width for each sample is different, which is very likely caused by an inhomogeneous sodium concentration and partial hydration that affect the superconducting state of $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_21.3\mathrm{H}_2\mathrm{O}`$. The $`\mathrm{T}_\mathrm{c}`$ values vary with Na content, showing the highest $`\mathrm{T}_\mathrm{c}`$ at 4.9K with *x*=0.42 given in the inset of Fig. 10. As shown in the inset the onset of the superconducting transition is identical with the divergence of FC and ZFC data.

Figure 11 is the plot of $`\mathrm{T}_\mathrm{c}`$ as function of sodium concentration of $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_21.3\mathrm{H}_2\mathrm{O}`$. For samples with *x*$``$0.22 and *x*$``$0.47 there is no superconducting transition detectable. These data agree to some extend with a recent study that shows a constant $`\mathrm{T}_\mathrm{c}`$ for Na contents up to 0.37 milne04 . However, they conflict with the superconducting ”dome” of $`\mathrm{T}_\mathrm{c}`$(x) demonstrated for powder samples 3 . This latter dependence motivated the proposal of intrinsic critical concentrations *x*cr1$``$1/4 and $`x_{\mathrm{cr2}}`$$``$1/3 related to charge ordering instabilities of the $`\mathrm{Co}^{3+}/\mathrm{Co}^{4+}`$ baskaran04 . In the present single crystal study the upper limit is shifted to much higher concentrations, *x*cr2 $``$ 0.45. This implies, that a charge ordering at *x*cr2$``$1/3 is not relevant to suppress $`\mathrm{T}_\mathrm{c}`$.

However, the superconducting phase extends now close to the phase line where non-hydrated $`\mathrm{Na}_{0.5}\mathrm{CoO}_2`$ shows a metal-insulator transition huang04 . This is evidence for an importance of electronic correlations at higher *x* in the hydrated system. It is interesting to note that the Néel transition $`\mathrm{T}_\mathrm{N}`$$``$20K is also only observed for higher Na concentrations *x*$`>`$0.75 in non-hydrated crystals, and that this ordering temperature does not change appreciably with *x*. Bulk antiferromagnetic order has been proven using myon spin rotation and other thermodynamic experimental technique 6 . The exact nature of the ordering, however, is still needed to further investigate.

Evidence for superlattice formation and electronic instabilities either due to Na or $`\mathrm{Co}^{3+}/\mathrm{Co}^{4+}`$ charge ordering are recently accumulating for non-hydrated $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ with *x*=0.5, but also to some extend for other stoichiometries zandbergen04 ; gavilano04 . Hydration contributes to the charge redistribution due to its effect on the $`\mathrm{CoO}_2`$ layer thickness and the formation of $`\mathrm{Na}(\mathrm{H}_2\mathrm{O})_4`$ clusters that shield disorder of the partially occupied Na sites. The above mentioned instabilities should influence the superconducting state as they modulate the electronic density of states and change nesting properties of the Fermi surface tanaka03 ; johannes04 . Further studies are needed to investigate a possible interrelation of Na content, hydration and structural details on well-characterized samples with highest $`\mathrm{T}_\mathrm{c}`$.

## IV Conclusion

Using TG-DTA we found that $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ is an incongruently melting compound, and single crystals can be grown using the optical floating zone technique. Both non-hydrated and fully hydrated crystals are exceptionally sensitive under ambient conditions. Study of XRD and TG indicates that the semi-hydrated phase $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_20.6\mathrm{H}_2\mathrm{O}`$ is more stable than the other compositions under ambient conditions. The XRD patterns show the presence of cells corresponding to the coexistence of several different hydrated phases in rehydrated samples. The superconducting transition temperature of $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_21.3\mathrm{H}_2\mathrm{O}`$ is only weakly influenced by the Na content. Highest $`\mathrm{T}_\mathrm{c}`$4.9 K is found for x$``$0.42, i.e. in the proximity of a metal-insulator transition observed in the non-hydrated phase with *x*=0.5. The hydrated compound is strongly unstable concerning its chemical, structural and electrical properties. We therefore propose to carefully reinvestigate the recently claimed evidence for unconventional superconductivity on single crystals.

Acknowledgements: We thank E. Brücher, G. Götz, E. Peters, and E. Winckler for technical support.


format: text

text:
# Core-collapse supernovae and gravitational waves

## 1 CORE-COLLAPSE SUPERNOVAE

Core-collapse supernovae—those of Type Ib, Ic, and II—result from the catastrophic collapse of the core of a massive star. Depending on the properties of the progenitor star, the collapse may result in a black hole; alternatively, the remnant can be a neutron star with $`GM/R0.1c^2`$, where $`G`$ is the gravitational constant, $`M`$ and $`R`$ are the neutron star mass and radius, and $`c`$ is the speed of light. Such a system is manifestly relativistic; and if the high densities and infall velocities implied by the above relation are combined with sufficient asphericity, the violence of core collapse and its aftermath may be expected to produce significant gravitational radiation. (Type Ia supernovae, which are caused by the thermonuclear detonation or deflagration of a white dwarf star, are not expected to be interesting sources of gravitational radiation.)

I begin by outlining our current understanding (and lack of understanding) of the core-collapse supernova process.

For most of their existence stars burn hydrogen into helium. In stars at least eight times as massive as the Sun, temperatures and densities become sufficiently high to burn to carbon, oxygen, neon, magnesium, and silicon and iron group elements. The iron group nuclei are the most tightly bound, and here burning in the core ceases.

The iron core—supported by electron degeneracy pressure—eventually becomes unstable. Its inner portion undergoes homologous collapse (velocity proportional to radius), and the outer portion collapses supersonically. Electron capture on nuclei is one instability leading to collapse, and this process continues throughout collapse, producing neutrinos. These neutrinos escape freely until densities and temperatures in the collapsing core become so high that even neutrinos are trapped.

Collapse is halted soon after the matter exceeds nuclear density; at this point (called “bounce”), a shock wave forms at the boundary between the homologous and supersonically collapsing regions. The shock begins to move out, but after the shock passes some distance beyond the surface of the newly born neutron star, it stalls as energy is lost to neutrino emission and dissociation of heavy nuclei falling through the shock.

The details of how the stalled shock is revived sufficiently to continue plowing through the outer layers of the progenitor star are unclear. Some combination of neutrino heating of material behind the shock, convection, instability of the spherical accretion shock, rotation, and magnetic fields launches the explosion.

It is natural to consider neutrino heating as a mechanism for shock revival, because neutrinos dominate the energetics of the post-bounce evolution. Initially, the nascent neutron star is a hot thermal bath of dense nuclear matter, electron/positron pairs, photons, and neutrinos, containing most of the gravitational potential energy released during core collapse. Neutrinos, having the weakest interactions, are the most efficient means of cooling; they diffuse outward on a time scale of seconds, and eventually escape with about 99% of the released gravitational energy.

Because neutrinos dominate the energetics of the system, a detailed understanding of their evolution will be integral to any detailed and definitive account of the supernova process. If we want to understand the explosion—which accounts for only about 1% of the energy budget of the system—we should carefully account for the neutrinos’ much larger contribution to the energy budget.

What sort of computation is needed to follow the neutrinos’ evolution? Deep inside the newly-born neutron star, the neutrinos and the fluid are tightly coupled (nearly in equilibrium); but as the neutrinos are transported from inside the neutron star, they go from a nearly isotropic diffusive regime to a strongly forward-peaked free-streaming region. Heating of material behind the shock occurs precisely in this transition region, and modeling this process accurately requires tracking both the energy and angle dependence of the neutrino distribution functions at every point in space.

While a full treatment of this six-dimensional neutrino radiation hydrodynamics problem remains too costly for currently available computational resources, there is much that has been learned over the years through detailed modeling.

## 2 SIMULATING THE EXPLOSION

Supernovae have a rich phenomenology—observations of many types that modelers would like to reproduce and explain. Chief among these is the explosion itself, which is not yet produced robustly and convincingly in simulations. Other observables of interest include neutrino signatures; neutron star spins, kick velocities, and magnetic fields; synthesized element abundances; all kinds of measurements across the electromagnetic spectrum; and of course the subject of this conference session, gravitational waves.

Simulations of core collapse, bounce, and its immediate aftermath have mostly aimed at the first few of these observables: the explosion mechanism, neutrino signatures, remnant pulsar properties, and gravitational waves. I will now describe some of the progress in this work in the past decade or so, focusing in particular on the explosion mechanism.

Throughout the 1990s, several groups performed simulations in two spatial dimensions. Even in two spatial dimensions, computational limitations required approximations that simplified the neutrino transport.

One simplification allowed for neutrino transport in two spatial dimensions, but with neutrino energy and angle dependence integrated out—effectively reducing a five dimensional problem to a two dimensional one (see for example ). These simulations exhibited explosions, suggesting that the enhancements in neutrino heating behind the shock resulting from convection provided a robust explosion mechanism. More recent simulations in three spatial dimensions with this same approximate treatment of neutrino transport showed similar outcomes .

A different simplification of neutrino transport employed in the 1990s was the imposition of energy-dependent neutrino distributions from spherically symmetric simulations onto fluid dynamics computations in two spatial dimensions . Unlike the multidimensional simulations discussed above, these did not exhibit explosions, casting doubt upon claims that convection-aided neutrino heating constituted a robust explosion mechanism.

The nagging qualitative difference between multidimensional simulations with different neutrino transport approximations renewed the motivation for simulations in which both the energy and angle dependence of the neutrino distributions were retained. Of necessity, the first such simulations were performed in spherical symmetry (actually a three-dimensional problem, depending on one space and two momentum space variables). Results from three different groups are in accord: Spherically symmetric models do not explode, even with solid neutrino transport .

Recently, one of these groups performed simulations in two spatial dimensions, in which their energy- and angle-dependent neutrino transport was made partially dependent on spatial polar angle as well as radius . Explosions were not seen in any of these simulations, except for one in which certain terms in the neutrino transport equation corresponding to Doppler shifts and angular aberration due to fluid motion were dropped. This was a surprising qualitative difference induced by terms contributing what are typically thought of as small corrections. The continuing lesson is that getting the details of the neutrino transport right makes a difference.

Where, then, do simulations aiming at the explosion mechanism stand? The above history suggests that elucidation of the mechanism will require simulations that feature truly spatially multidimensional neutrino transport. In addition, inclusion of magnetic field dynamics seems increasingly strongly motivated as a possible driver of the explosion, because simulations with “better” neutrino transport have failed to explode—even in multiple spatial dimensions.

## 3 GRAVITATIONAL RADIATION

Core bounce and associated phenomena (halt of the collapse, shock formation, and “ring-down”) would be the strongest source of gravitational radiation from a core-collapse supernova. Should such an event occur in our galaxy (this happens only once every $``$few-several decades), it probably would be detectable by first-generation gravitational wave interferometers. This conclusion is reached in several studies; a recent example is the work of Ott et al. .

These authors also make some points that highlight the need for more complete simulations. Most often, the results of spherically symmetric stellar evolution computations are artificially put into rotation for use as progenitors in studies of core collapse; but in simulations with the most recent progenitors that include rotation and magnetic fields, the gravitational radiation can be an order of magnitude weaker. In addition, some features of the waveform (such as coherent post-bounce oscillations) probably will change with the inclusion of realistic neutrino transport because of the associated stalling of the shock.

In addition to the violence of core bounce and its immediate aftermath, additional phenomena in the subsequent evolution may produce gravitational radiation at lower amplitudes. These include triaxial instabilities in the neutron star, anisotropic neutrino emission, convection, and (on longer time scales than the supernova process itself) R-mode instabilities . Another less well-known possible source of gravitational waves is the standing accretion shock instability .


text:
# Effect of Differential Rotation on the Maximum Mass of Neutron Stars: Realistic Nuclear Equations of State

## 1 Introduction

The remnant formed in the coalescence of binary neutron stars is likely to be differentially rotating (see, e.g., the dynamical simulations of Rasio & Shapiro 1992, 1994, 1999; Shibata & Uryū 2000, 2002; Faber, Rasio & Manor 2001; Oechlin, Rosswog & Thielemann, 2002; Shibata, Taniguchi & Uryū 2003; Faber, Grandclément & Rasio 2003; see also the review of Baumgarte and Shapiro 2003). It is likely that differential rotation will play an important role in the dynamical stability of these remnants, since it can be very effective in increasing the their maximum allowed mass (Baumgarte, Shapiro & Shibata 2000, hereafter BSS; Lyford, Baumgarte & Shapiro 2003, hereafter LBS).

Most neutron stars in binaries have individual gravitational masses close to 1.4 $`M_{}`$ (compare Table 2 below). Furthermore, most recent realistic nuclear equations of state predict a maximum allowed mass for nonrotating neutron stars in the range of about 1.7 – 2.3 $`M_{}`$ (compare Table 1 below). Taken together, these facts seem to suggest that the coalescence of binary neutron stars would lead to prompt collapse to a black hole. However, rotation, and especially differential rotation, can increase the maximum allowed mass significantly111For coalescence from the innermost stable circular orbit, shear and shock-induced thermal pressure, which can also increase the maximum allowed mass, is found to have a smaller effect..

The maximum allowed mass of uniformly rotating stars is limited by the spin rate at which the fluid at the equator moves on a geodesic (the Kepler limit); any further speed-up would lead to mass shedding. Uniform rotation can therefore increase the maximum allowed mass by about 20 % at most for very stiff equations of state (Cook, Shapiro & Teukolsky, 1992, 1994a, 1994b, hereafter CST1, CST2 and CST3 respectively; see also Table 1 below), which is not sufficient to stabilize remnants of binary neutron star coalescence. Uniformly rotating neutron stars with rest masses exceeding the maximum allowed rest mass for non-rotating stars (for the same equation of state) are referred to as supramassive neutron stars.

Differential rotation, however, is much more effective in increasing the maximum allowed mass. Unlike a uniformly rotating star, the rotation rate at the core of a differentially rotating star is not restricted to the maximum rotation rate at the equator, so that the core can be supported by rapid rotation without the equator having to exceed the Kepler limit. This effect was demonstrated in Newtonian gravitation by Ostriker, Bodenheimer & Lynden-Bell (1966) for white dwarfs, and in general relativity by BSS for $`n=1`$ polytropes. BSS also showed by way of illustration that stars with about 60% more mass than the maximum allowed mass of the corresonding nonrotating star can be dynamically stable against both radial and nonaxisymmetric modes. BSS refer to differentially rotating equilibrium configurations with rest masses exceeding the maximum rest mass of a uniformly rotating star as hypermassive neutron stars. LBS generalized the equilibrium results of BSS to other polytropic indices in the range $`0.5n2.9`$, and found that the largest relative increases in the maximum allowed mass can be found for polytropic indices close to $`n=1`$. It is therefore not surprising that the hypermassive binary neutron star remnants formed in the $`n=1`$ simulations of Shibata, Taniguchi & Uryū (2003) do not collpase promptly to black holes unless they exceed the maximum nonrotating mass by more than about 70%. The merged remnants will be dynamically stable on a dynamical timescale. Viscous damping and magnetic braking of differential rotation will likely occur on a secular timescale (which is much greater than the dynamical timescale), leading to a delayed collapse and a delayed burst of gravitational radiation (BSS; Shapiro 2000; Cook, Shapiro & Stevens 2003; Liu & Shapiro 2003).

In this paper, we generalize the results of BSS and LBS to realistic equations of state. We construct differentially rotating neutron stars for a sample of six cold nuclear equations of state, and tabulate their maximum allowed masses as a function of the degree of differential rotation. In addition, we identify an “effective” polytropic index for each of these equations of state, and construct differentially rotating polytropes for these indices. We consistently find larger mass increases for the polytropes than for the nuclear equations of state, which is due to a drop in the maximum density for a rotating star with a large degree of differential rotation. However, even for nuclear equations of state we find increases in the maximum allowed rest mass on the order of about 50%, so that binary neutron star coalescence should often result in hypermassive neutron stars. We explicitly predict the fate of the six double neutron star binaries that have been identified to date, using the results of our equilibrium model calculations.

The paper is organized as follows. In Section 2 we briefly summarize the numerical method and discuss the nuclear equations of state adopted in our calculations. In Section 3 we present numerical results and compare our findings for nuclear equations of state with those for polytropes. We summarize our findings and discuss consequences for the six observed binary neutron star systems in Section 4. In Appendix A we tabulate our numerical results. We adopt gravitational units and set $`G=c=1`$.

## 2 Constructing Numerical Models

### 2.1 Equilibrium Models of Differentially Rotating Stars

As in BSS and LBS, we use a modified version of the numerical code of CST to construct equilibrium models of differentially rotating neutron stars in general relativity. The code is based on similar numerical methods developed by Hachisu (1986) and Komatsu, Eriguchi, & Hachisu (1989), and we refer to CST1 for details.

Constructing differentially rotating neutron star models requires choosing a rotation law $`F(\mathrm{\Omega })=u^tu_\varphi `$, where $`u^t`$ and $`u_\varphi `$ are components of the four-velocity $`u^\alpha `$ and $`\mathrm{\Omega }`$ is the angular velocity. We follow CST1 and assume a rotation law $`F(\mathrm{\Omega })=A^2(\mathrm{\Omega }_c\mathrm{\Omega })`$, where the parameter $`A`$ has units of length and where $`\mathrm{\Omega }_c`$ is the central value of the angular velocity. Expressing $`u^t`$ and $`u_\varphi `$ in terms of $`\mathrm{\Omega }`$ and metric potentials yields equation (42) of CST1, or

$$\mathrm{\Omega }=\frac{\mathrm{\Omega }_c}{1+\widehat{A}^2\widehat{r}^2\mathrm{sin}^2\theta }$$
(1)
in the Newtonian limit. Here we have rescaled $`A`$ and $`r`$ in terms of the equatorial radius $`R_e`$: $`\widehat{A}A/R_e`$ and $`\widehat{r}r/R_e`$. The parameter $`\widehat{A}`$ is a measure of the degree of differential rotation and determines the length scale over which $`\mathrm{\Omega }`$ changes. Since uniform rotation is recovered in the limit $`\widehat{A}^10`$, it is convenient to parametrize sequences by $`\widehat{A}^1`$. In the Newtonian limit the ratio between the central and equatorial angular velocities $`\mathrm{\Omega }_c/\mathrm{\Omega }_e`$ is related to $`\widehat{A}`$ by $`\mathrm{\Omega }_c/\mathrm{\Omega }_e=1+\widehat{A}^2`$, but for relativistic configurations this relation holds only approximately.

We adopt this particular rotation law for convenience and for easy comparison with many other authors who have assumed the same law. However, as we pointed out in LBS, this rotation law approximates reasonably well the angular velocity profile of binary neutron star remnants in the fully relativistic dynamical simulations of Shibata & Uryu (2000, 2002) and the post-Newtonian simulations of Faber, Rasio, & Manor (2001; see also Faber, Grandclément & Rasio 2003). This suggests that the above rotation law may provide a reasonable parametrization of differential rotation profiles that one might expect to find in binary neutron star merger remnants.

We modified the numerical algorithm of CST1 by fixing the maximum interior density instead of the central density for each model. This change allows us to construct higher mass models in some cases, since the central density does not always coincide with the maximum density and hence may not specify a model uniquely.

For each equation of state and a given value of $`\widehat{A}`$ we construct a sequence of models for each value of the maximum density by starting with a static, spherically symmetric star and then decreasing the ratio of the polar to equatorial radius, $`R_p/R_e`$, in decrements of 0.015. This sequence ends when we reach mass shedding (for large values of $`\widehat{A}`$) or when the code fails to converge (indicating the termination of equilibrium solutions) or when $`=0`$ (beyond which the star would become a toroid). For each one of these sequences, the maximum achieved mass is recorded. We repeat this procedure for different values of the maximum density, which yields the maximum achieved mass as a function of maximum density. The maximum of this curve is the maximum allowed mass for this particular equation of state and the chosen value of $`\widehat{A}`$. In the Appendix we tabulate our numerical results for the equations of state described in Section 2.2, and for $`\widehat{A}^1=0,0.3,0.5,0.7`$ and 1.0. Clearly, our maximum allowed masses are lower limits in the sense that even higher mass models may exist for other values of $`\widehat{A}`$ or different differential rotation laws.

### 2.2 Equations of State

In this paper we adopt six cold, nuclear equations of state (EOSs) which are listed, together with numerical results, in Table 1. The first five, A, D, L, UT and FPS are adopted from CST3 (who in turn adopted the labeling of Friedman, Ipser, & Parker (1986) for the first three EOSs). The last EOS, APR, was adapted from Akmal, Pandharipande & Ravenhall (1998) by Ravenhall and Pandharipande (with the smooth matching between low and high density provided by G. Cook). In Figure 1 we plot the pressure $`P`$ as a function of the rest-mass density $`\rho _0`$ for these EOSs.

EOS A (Pandharipande 1971) models the interaction of neutrons at high densities with a Reid soft-core potential. EOS D is model V of Bethe & Johnson (1974). In EOS L the nucleon interaction is modeled in terms of a mean scalar field (Pandharipande, & Smith 1975). Both EOSs UT (Wiringa, Fiks, & Fabrocini 1988) and FPS (Lorenz, Ravenhall, & Pethick, 1993) are modern versions of an earlier EOS proposed by Friedman & Pandharipande (1981), which employs both two-body (U14) and three-body nucleon interactions (TNI). EOS UT improves the treatment of matter at high densities, while FPS describes the interactions in terms of a Skyrme model. EOS APR (Akmal, Pandharipande, & Ravenhall 1998) adopts a modern two-nucleon interaction (A18) together with boost corrections, as well as the UIX three-nucleon potential.

The different descriptions of nucleon interactions affect the EOS only at high densities. At low densities ($`\rho _010^4\text{g cm}^3`$), the EOSs employ the Feynman, Metropolis, & Teller (1949) EOS, joining onto the Baym, Pethick, & Sutherland (1971) EOS up to neutron drip at $`\rho _04\times 10^{11}\text{g cm}^3`$. Above neutron drip, EOSs A, D join onto the Baym, Bethe, & Pethick (1971) EOS, while EOSs L and UT join onto the Negele & Vautherin (1973) EOS. EOS APR joins onto the FPS EOS below a number density of 0.1 fm-3 ($`\rho _0=1.26\times 10^{14}\text{g cm}^3`$),

All EOSs are read into our numerical code in tabular form, listing rest-density $`\rho _0`$, the total energy density $`ϵ`$, and the pressure $`P`$ at discreet points. Intermediate values of these quantities are computed by interpolation. Further details on the numerical implementation of these EOSs, including the fitting and tabulation, can be found in CST3.

## 3 Results

Our numerical results are summarized in Table 1, where we list, for each equation of state, the nonrotating (Tolman-Oppenheimer-Volkoff, TOV) maximum masses (both gravitational mass $`M^{\mathrm{TOV}}`$ and rest mass $`M_0^{\mathrm{TOV}}`$), the maximum masses for uniform rotation, and the maximum masses that we found in our survey of differential rotation. For convenience we also list fractional mass increases with respect to the TOV maximum mass, as well as the parameter $`\widehat{A}^1`$ for which the maximum mass was encountered. More detailed results, with results for all values of $`\widehat{A}^1`$ considered in this paper, are tabulated in Appendix A.

The fractional mass increases are typically in the order of 50%. These increases are still large enough to support the remnant of a coalescing binary neutron star system in many cases (compare Section 4). However, they are noticeably below the increases for polytropic EOSs. LBS found that the fractional mass increases depend strongly on the polytropic index $`n`$, i.e. the stiffness of the EOS. For values of $`n`$ between 0.75 and 1.25, corresponding to moderately stiff EOSs, they found relative mass increases exceeding 100%. It is therefore somewhat surprising that the nuclear EOSs considered in this paper yield significantly smaller numbers.

To understand this result we construct models of differentially rotating neutron stars for polytropic EOSs

$$P=K\rho _0^{1+1/n},$$
(2)
where $`K`$ a polytropic constant and $`n`$ the polytropic index. For each nuclear EOS we identify an “effective” polytropic index $`n_{\mathrm{eff}}`$ as follows. We first construct, for a nuclear EOS, the nonrotating (TOV) maximum-mass model, and compute for this model the central concentration $`(ϵ_c/\overline{ϵ})_{\mathrm{nucl}}`$. Here $`ϵ_c`$ is the central value of the energy density $`ϵ`$, and we define the average energy density from

$$\overline{ϵ}\frac{3M}{4\pi R^3},$$
(3)
where $`M`$ is the gravitational mass and $`R`$ the circumferential radius. In Newtonian gravitation, this central condensation would correspond to a unique value of $`n`$. In general relativity, however, this is no longer the case. To identify the effective polytropic index $`n_{\mathrm{eff}}`$, we therefore construct TOV maximum mass configurations for a large sample of polytropic indices $`n`$, and compute their central condensation $`(ϵ_c/\overline{ϵ})_{\mathrm{poly}}`$ as described above (compare Table 1 in LBS). We then interpolate between the polytropic values of the central condensation $`(ϵ_c/\overline{ϵ})_{\mathrm{poly}}`$ to the desired nuclear value $`(ϵ_c/\overline{ϵ})_{\mathrm{nucl}}`$, and thereby identify the effective polytropic index $`n_{\mathrm{eff}}`$. We fix the polytropic constant $`K`$ in (2) by setting the mass of the polytropic maximum mass TOV model, which scales with $`K^{n/2}`$, equal to its counterpart for the corresponding nuclear EOS.

We consistently find that the agreement between the nuclear EOS models and the corresponding polytropes is quite good (to within 10% or so) for those models that have a maximum density $`ϵ_{\mathrm{max}}`$ similar to the maximum mass TOV star. However, for larger degrees of differential rotation the maximum mass models shift to lower values of $`ϵ_{\mathrm{max}}`$ as the star become increasingly “toroidal” and bloated perpendicular to the rotation axis. We show this effect for EOS FPS in Figure 2, and we find qualitatively indentical results for all other EOSs (see also Figure 1 in BSS).

This behavior can be understood by computing, as an alternative to our construction above, a “nuclear” polytropic index from the slope of the graphs in Fig. 1. Using the adiabatic index

$$\mathrm{\Gamma }=\frac{\mathrm{ln}P}{\mathrm{ln}\rho _0}$$
(4)
we difine a nuclear polytropic index by

$$n_{\mathrm{nucl}}=\frac{1}{\mathrm{\Gamma }1}.$$
(5)
Clearly this index is a function of density, and its value reflects different interactions dominating different density regimes. Across phase transitions, $`n_{\mathrm{nucl}}`$ can even change discontinously. We include a plot of $`n_{\mathrm{FPS}}`$ as a function of the density $`ϵ`$ in the inserted panel in Figure 2.

For a polytropic EOS, on the other hand, $`n`$ is strictly independent of the density. We can therefore expect good agreement between polytropic and nuclear models only as long as $`n_{\mathrm{eff}}`$ is close to $`n_{\mathrm{nucl}}`$ for most of the matter in the star. As can be seen in the panel in Figure 2, $`n_{\mathrm{eff}}`$ is very close to $`n_{\mathrm{FPS}}`$ for densities in the range $`(1.53)\times 10^{15}\text{g}\text{cm}^3`$. This coincides quite well with the regime of maximum densities in equilibrium models for which the agreement between the polytropic and nuclear models is fairly good. For lower densities, $`n_{\mathrm{FPS}}`$ drops significantly below $`n_{\mathrm{eff}}`$, and it is therefore not surprising that stellar models with maximum densities in this regime no longer agree very well.

For all other EOSs we find very similar results. In all cases, the nuclear adiabatic index $`n_{\mathrm{nucl}}`$ is very close to the effective polytropic index $`n_{\mathrm{eff}}`$ (as determined from the central condensation) for densities close the the maximum density of the maximum mass TOV star. Matter in this density regime dominates the structure of the star for uniform and moderate degrees of differential rotation, leading to good agreement between nuclear and polytropic models. For larger degrees of differential rotation, the star becomes increasingly “toroidal” and bloated, and the maximum density shifts to smaller values. As a consequence, the star is dominated by matter at smaller densities, for which $`n_{\mathrm{nucl}}`$ no longer agrees with $`n_{\mathrm{eff}}`$. Accordingly, the results for the polytropic models are no longer in good agreement with the nuclear models.

## 4 Summary and Discussion

We have constructed models of differentially rotating neutron stars in general relativity for a sample of realistic nuclear equations of state. We find that their maximum rest masses are typically larger than those of non-rotating neutron stars by about 50%. This increase is significantly less than the increases for moderately stiff polytropic equations of state (compare LBS). This deviation can be explained in terms of a drop in the maximum density for stars with large degrees of differential rotation. At these smaller densities, the stiffness of the nuclear equation of state, as determined from equation (5), is different from that at the maximum density of maximum mass nonrotating or uniformly rotating stars. Consequently, predictions based on simple polytropic models break down for sufficiently large degrees of differential rotation.

Even a mass increase of 50% may well be sufficient to support the remnant of binary neutron star coalescence. In fact, our findings can be used to predict possible fates for the six known double neutron star binaries as listed in Table 2. For four of those, B1524+12, B1913+16, B2127+11C and the newly discovered binary J0737-3039, the individual gravitational masses are known to high accuracy (Thorsett & Chakrabarty 1999; Burgay et. al. 2003, Lyne et. al. 2004). For the two remaining ones, J1518+4904 and B2303+46, the orbital inclination is unknown, so that the combined gravitational mass is known to much higher accuracy than the individual masses. Assuming no mass loss during the coalescence and merger (Faber, Grandclément & Rasio 2003) the total rest mass of the binary is conserved (the total gravitational mass is not conserved, due to neutrino and gravitational radiation losses). The rest mass corresponding to a particular gravitational mass depends on the equation of state, as is shown in Figure 3 for the four binaries with well-established individual masses.

The total rest mass $`M_0^{\mathrm{tot}}`$ of the binary is an indicator of the fate of the merger remnant. If the total rest mass is less than the maximum (nonrotating) TOV mass $`M_0^{\mathrm{TOV}}`$, then the remnant can be supported against gravitational collapse even without rotation. This possibility is identified by “NS” in Table 2. If $`M_0^{\mathrm{tot}}`$ is greater than $`M_0^{\mathrm{TOV}}`$, but less than the maximum mass of uniformly rotating stars $`M_0^{UNI}`$, then the merger could lead to a supramassive star, supported by uniform rotation (denoted by “SNS” in Table 2). If $`M_0^{\mathrm{tot}}`$ is greater than $`M_0^{\mathrm{UNI}}`$ but less than the maximum mass we find for differentially rotating stars $`M_0^{\mathrm{DIF}}`$, merger may lead to a hypermassive neutron star supported by differential rotation (denoted by “HNS”). If, finally, $`M_0^{\mathrm{tot}}`$ is greater than $`M_0^{\mathrm{DIF}}`$, the remnant will promply collapse to a black hole (“BH”). Clearly, these limits are only approximations, and the true outcome will depend not only on the mass and EOS but also on the non-equilibrium dynamics of the merger.

Since the rest masses of the observed binaries depend on the equation of state, so do the possible fates. Moreover, since all binary neutron star systems known to date have similar total masses, their fate, for a given equation of state, is presumably similar. In Table 2, the predictions for four of the six binaries are identical; only the high mass binary B1913+16 and the low mass binary J0737-3039 show any differences for some of the equations of state.

Table 2 suggests two immediate conclusions. For all the recent equations of state we survey and for all binaries except for B1913+16 the coalescence is likely to lead to a neutron star. For three out of the six equations of state, the merger of these binaries could lead to a hypermassive neutron star. Unlike supramassive neutron stars, hypermassive neutron stars last only a fairly short time, until magnetic fields or viscous dissipation drive the objects toward uniform rotation, triggering a delayed collapse to a black hole (see BSS).

Secondly, Table 2 shows that the outcome of binary neutron star coalescence is quite sensitive to the equation of state. This means that the observational identification of a merger remnant, and its possible delayed collapse to a black hole, could place a tight constraint on the equation of state. Advanced LIGO (the Laser Interferometer Gravitational Wave Observatory) and other new gravitational wave laser interferometers may be able to observe binary neutron star coalescence. Comparing the most recent equations of state, UT, FPS and APR, it seems likely that the coalescence will lead to a rotationally supported neutron star. Observation of the remnant’s delayed collapse to a black hole, or its absence, would then help to put contraints on possible equations of state.

It is a pleasure to thank G. Ravenhall and V. J. Pandharipande for providing the data for equation of state APR. This paper was supported in part by NSF Grant PHY-0139907 at Bowdoin College and NSF Grants PHY-0090310 and PHY-0205155 and NASA Grant NAG 5-10781 at the University of Illinois at Urbana-Champaign.

## Appendix A Numerical Results for Maximum Masses

We list in the Tables below values for the maximum rest masses for uniformly and differentially rotating polytropes. For each equation of state we tabulate the differential rotation parameter $`\widehat{A}^1`$, the maximum energy density $`ϵ_{\mathrm{max}}`$, the ratio of the central and equatorial angular velocity $`\mathrm{\Omega }_c/\mathrm{\Omega }_e`$, the ratio between polar and equatorial radii $`R_p/R_e`$, the ratio of the (relativistic) rotational kinetic energy and the gravitational binding energy $`T/|W|`$, the maximum rest mass $`M_0^{\mathrm{max}}`$, and the fractional rest mass increase $`(\delta M_0/M_0)`$. We also construct differentially rotating polytropes for effective polytropic indices $`n_{\mathrm{eff}}`$ as described in Section 3 and listed in each Table heading. We include in the Tables the fractional mass increases $`(\delta M_0/M_0)_{\mathrm{poly}}`$ for these polytropes.

All models are computed with 64 zones both in the radial and angular direction, and with the Legendre polynomial expansion truncated at $`\mathrm{}=16`$ (see CST1 for details of the numerical implementation). The error in the ratio $`R_p/R_e`$ is determined by our stepsize of 0.015 in this ratio (compare Section 2.1). We finally note that highly toroidal configurations depend very sensitively on the input parameters, so that those mass increases depend on the parameters adopted in our restricted survey.


text:
# Quantum theory of shuttling instability in a movable quantum dot array

## 1 Introduction

Since the first proposal by Gorelik et al.  of the shuttling instability in a generic nanoelectromechanical system (NEMS) consisting in a movable single electron transistor the phenomenon of the shuttling transport attracted much attention. However, until recently the fully quantum theory of the phenomenon was not developed. The first quantum mechanical study on a modified setup was accomplished by Armour and MacKinnon , closely followed by the work of the present authors on the original Gorelik’s setup . In this paper we extend the results by Armour and MacKinnon. The phase space analysis in terms of the Wigner functions introduced in  reveals directly the nature of the transport (incoherent tunnelling versus shuttling) in different regions in contrast to the indirect evidence used by Armour and MacKinnon. Also the new numerical scheme that we use enables us to access a wider range of parameters fast and reliably.

## 2 Model and method of solution

Let us consider a simple NEMS consisting of an array of three quantum dots (device) connected to two leads . The central dot is assumed to be movable in a parabolic potential. The Hamiltonian of the device consists of the mechanical and the electronic parts $`H=H_{\mathrm{mech}}+H_{\mathrm{el}}`$ where $`H_{\mathrm{mech}}=\mathrm{}\omega \widehat{a}^{}\widehat{a}`$ is the Hamiltonian of the harmonic oscillator and $`H_{\mathrm{el}}=_{\alpha ,\beta \{0,L,C,R\}}|\alpha ϵ_{\alpha \beta }(\widehat{x})\beta |`$. We assume strong Coulomb blockade regime with no double occupancies so that the vectors $`|\alpha `$ with $`\alpha =0,L,C,R`$ span the entire electronic Hilbert space of the device. Each matrix element $`ϵ_{\alpha \beta }(\widehat{x})`$ is still a full matrix in the oscillator space. $`ϵ(\widehat{x})`$ reads explicitly

$$ϵ(\widehat{x})=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& \frac{\mathrm{\Delta }V}{2}& t_L(\widehat{x})& 0\\ 0& t_L(\widehat{x})& \frac{\mathrm{\Delta }V}{2x_0}\widehat{x}& t_R(\widehat{x})\\ 0& 0& t_R(\widehat{x})& \frac{\mathrm{\Delta }V}{2}\end{array}\right]$$
(1)
where $`\widehat{x}=\sqrt{\frac{\mathrm{}}{2m\omega }}(\widehat{a}+\widehat{a}^{})`$ is the position operator, $`\mathrm{\Delta }V`$, the device bias, is the difference between the energy of the left and the right dot. $`x_0`$ is half the distance between the two outer dots and represents the maximum amplitude of the inner dot oscillation. The three dots are electrically connected only via a tunnelling mechanism. The tunnelling length is given by $`1/\alpha `$ and the tunnelling strengths depend on the position operator $`\widehat{x}`$ of the inner grain as $`t_L(\widehat{x})=V_0e^{\alpha (x_0+\widehat{x})},t_R(\widehat{x})=V_0e^{\alpha (x_0\widehat{x})}`$.

The dynamics of the device is described by the generalized master equation (GME) formalism . The density matrix of the device evolves according to $`\dot{\rho }=i[H,\rho ]+\mathrm{\Xi }\rho +\dot{\rho }_d`$. The first term represents the coherent evolution of the isolated device. The coupling to the leads responsible for the electronic transfer from/to the device is introduced in the wide band approximation following Gurvitz et al.  and yields the second term in the equation (in the block notation of )

$$\mathrm{\Xi }\rho =\mathrm{\Gamma }\left[\begin{array}{cccc}\rho _{RR}\rho _{00}& 0& 0& 0\\ 0& \rho _{00}& 0& \rho _{LR}/2\\ 0& 0& 0& \rho _{CR}/2\\ 0& \rho _{RL}/2& \rho _{RC}/2& \rho _{RR}\end{array}\right]$$
(2)
with $`\mathrm{\Gamma }`$ being the injection rate to the leads. The third term describes the effect of the environment on the oscillator, consisting in mechanical damping and random quantum and thermal excitation (Langevin force). It reads

$$\dot{\rho }_d=\frac{\gamma }{2}\overline{n}(aa^{}\rho 2a^{}\rho a+\rho aa^{})\frac{\gamma }{2}(\overline{n}+1)(a^{}a\rho 2a\rho a^{}+\rho a^{}a)$$
(3)
where $`\gamma `$ is the damping rate and $`\overline{n}`$ is the mean occupation number of the oscillator at temperature $`T`$.

The stationary version of the above GME was solved numerically after the oscillator Hilbert space was truncated at sufficiently large $`N`$ so that all dynamically excited states were contained within the basis. Utilizing the decoupling properties of the GME in the block notation resulted in the problem of finding a unique null space of a (super)matrix with the linear dimension $`10N^2`$ (with up to $`N=40`$). The Arnoldi iteration  with preconditioning was found to be superior to direct methods, such as singular value decomposition or inverse iteration. The preconditioning involving the inversion of the Sylvester part of the problem (a fast procedure) is crucial for the convergence of the method which is very fast and low memory consuming (the whole (super)matrix does not have to be stored in the memory).

## 3 Results and discussion

The resulting stationary density matrix was used to evaluate the mean value of current $`I=\mathrm{\Gamma }\mathrm{Tr}_{\mathrm{osc}}\rho _{RR}`$  and the phase space distribution of the charged central dot using the transformation into Wigner coordinates . These quantities as a function of the device bias $`\mathrm{\Delta }V`$ and for three different values of the injection rate $`\mathrm{\Gamma }`$ are plotted in the figures 1 and 2, respectively. It was found in  that the triple dot system exhibits different regimes of transport at different device biases. The current peaks at $`\mathrm{\Delta }Vn\omega `$ (see figure 1) were identified as effects of electromechanical resonances within the device. Yet, the different peaks may correspond to different physical mechanisms — while the peak around $`\mathrm{\Delta }V\omega `$ is mainly due to the incoherent oscillator-assisted tunnelling the peak at $`\mathrm{\Delta }V2\omega `$ reveals a clear shuttling component. This finding by Armour and MacKinnon based on indirect evidence of parametric dependencies of the current curves (e.g. the dependence of the current curve on the tunneling length $`1/\alpha `$) is confirmed by the direct inspection of the phase space distributions (see the first row of figure 2). The half-moon-like shape characteristic for shuttling transport  is only present around $`\mathrm{\Delta }V2\omega `$ while all other plots show just the fuzzy spot indicative of incoherent tunnelling. However, our direct criterion for detecting the shuttling regime reveals a close similarity between the resonances. For increasing injection rate we can see that the shuttling regime gradually sets in also in the vicinity of the first resonance peak. This reveals the crossover character of the onset of the shuttling instability found also previously . The sharp transition into the shuttling regime reported in semiclassical studies is smeared into the crossover due to the noise present in the system and properly accounted for in our approach.

The authors thank T. Eirola for indispensable advice concerning the numerical methods. Support of the grant 202/01/D099 of the Czech grant agency for one of us (T.N.) is also gratefully acknowledged.


text:
# A novel formulation of nonlocal electrostatics

## Abstract

The accurate modeling of the dielectric properties of water is crucial for many applications in physics, computational chemistry and molecular biology. This becomes possible in the framework of nonlocal electrostatics, for which we propose a novel formulation allowing for numerical solutions for the nontrivial molecular geometries arising in the applications mentioned before. Our approach is based on the introduction of a secondary field, $`\psi `$, which acts as the potential for the rotation free part of the dielectric displacement field $`𝐃`$. For many relevant models, the dielectric function of the medium can be expressed as the Green’s function of a local differential operator. In this case, the resulting coupled Poisson (-Boltzmann) equations for $`\psi `$ and the electrostatic potential $`\varphi `$ reduce to a system of coupled PDEs. The approach is illustrated by its application to simple geometries.

The theory of continuum electrostatics plays a major role in the determination of solvation free energies of atoms, ions, and biomolecules honig95 . Recent progress in its applicability to biological systems has been impressive: the electrostatic potentials of large biomolecules such as, e.g., microtubuli and ribosomes, can be determined baker01 . Unfortunately, the standard continuum approach ultimately becomes inaccurate when used to determine electrostatic properties on atomic scales simonson01 , as it is featureless, i.e., the correlation between solvent arrangements and the geometrical structure of biomolecular assemblies is not taken into account. On the other hand, continuum electrostatics is still much more efficient from a computational point of view than microscopic simulations based on, e.g., molecular dynamics (MD). Therefore, interest has risen recently in extensions of the theory of continuum electrostatics that allow to account for spatial variations of the dielectric behaviour of the solvent, in particular near boundaries hansen ; honig95 . Part of the motivation for such approaches stems from the field of protein docking, where a realistic and efficient modelling of solvent properties is essential docking .

Within the continuum theory of electrodynamics, spatial dispersion effects can be taken into account in an approach called ‘nonlocal electrostatics’ korn1 ; korn2 ; korn3 ; korn4 . It rests on the assumption of a linear relationship between the dielectric displacement field and the electric field mediated by a permittivity kernel depending on two spatial arguments,

$$𝐃(𝐫)=\epsilon _0𝑑𝐫^{}\epsilon (𝐫,𝐫^{})𝐄(𝐫^{})$$
(1)
where $`\epsilon (𝐫,𝐫^{})`$ is the dielectric permittivity tensor. Equivalently to eq.(1) one can express the nonlocal relationship in terms of the polarization fields korn4 .

While the theory of nonlocal electrostatics remains firmly embedded within the well-understood framework of Maxwell’s theory, it introduces a new characteristic length scale absent in local electrostatics: the correlation length $`\lambda `$ of the polarization correlations between the solvent molecules. This length sets the relevant scale for the deviation of the dielectric properties of the solvent from its bulk value. Thus, nonlocal electrostatics is a serious candidate for a more realistic description of solvent properties, provided it is also computationally tractable. It is here where the difficulties arise, however. The theory of nonlocal electrostatics, discussed in detail below, is technically considerably more demanding than local electrostatics, since it is usually formulated as a system of coupled integro-differential equations. Consequently, it has so far only been applied to idealized situations, and even then typically after introducing additional approximations in order to obtain analytical results levadny ; cherep . For complex geometries, the solution of the equations by numerical methods becomes a formidable task.

Thus, a reformulation of the equations of nonlocal electrostatics is needed in order to make the theory applicable to real-world problems. Here we present a scheme which allows to rewrite the theory in terms of a system of partial differential equations for the local fields which consequently makes it amenable to standard methods of numerical analysis. This derivation relies on two assumptions which are typically valid for e.g. the discussion of solvation problems for biomolecules: (i) the linearity of the relationship between dielectric response and the electric field, and (ii) the representation of the dielectric function in terms of Green functions of known differential operators. We first derive the set of equations, and then illustrate how to solve them on (simple) examples.

In the following we consider the situation given in Figure 1. A domain $`\mathrm{\Omega }`$ (which represents a molecule) is embedded in the solvent $`\mathrm{\Sigma }`$ which takes up the whole space except $`\mathrm{\Omega }`$. The surface of the embedded domain is denoted by $`\mathrm{\Gamma }`$. Within $`\mathrm{\Omega }`$, the dielectric properties are assumed to be local, i.e., $`\epsilon \epsilon _\mathrm{\Omega }`$, while in the solvent space $`\mathrm{\Sigma }`$, eq.(1) holds. On large scales, i.e. when $`|𝐫𝐫^{}|\lambda `$, the dielectric response in the solvent is again local with the bulk dielectric constant $`\epsilon _\mathrm{\Sigma }`$ (with, e.g., $`\epsilon _\mathrm{\Sigma }80`$ for water). In this situation the standard equations of nonlocal electrostatics read as

$$\mathrm{\Delta }\varphi _\mathrm{\Omega }=\frac{1}{\epsilon _\mathrm{\Omega }\epsilon _0}\varrho ,𝐫ϵ\mathrm{\Omega }$$
(2)
$$\epsilon _0_\mathrm{\Sigma }𝑑𝐫^{}\epsilon (𝐫,𝐫^{})^{}\varphi _\mathrm{\Sigma }(𝐫^{})=0,𝐫,𝐫^{}ϵ\mathrm{\Sigma }$$
(3)
where $`\varrho `$ is the density of fixed charges which are assumed to lie confined within $`\mathrm{\Omega }\backslash \mathrm{\Gamma }`$, which is usually the case for biomolecules in solution. This is no restriction on the validity of our approach: the existence of surface charges on $`\mathrm{\Gamma }`$ only slightly modifies the boundary conditions. Eqs.(2,3) are the Poisson equations for the geometry of Fig. 1. The primed symbol $`^{}`$ denotes the differentiation with respect to $`𝐫^{}`$. Without surface charges on $`\mathrm{\Gamma }`$, the boundary conditions to eqs.(2,3) for the normal (n) and tangential (t) components of the electric and dielectric displacement field on the boundary $`\mathrm{\Gamma }`$ are given by

$$𝐃_{\mathrm{\Omega },n}=𝐃_{\mathrm{\Sigma },n},𝐄_{\mathrm{\Omega },t}=𝐄_{\mathrm{\Sigma },t}$$
(4)
where, by virtue of eq.(3), the boundary condition for $`𝐃`$ is also nonlocal.

In order to step over from the integro-differential to a purely differential formulation we introduce, in addition to the potential field $`\varphi `$, another potential field $`\psi `$ within both compartments $`\mathrm{\Omega }`$ and $`\mathrm{\Sigma }`$. Attempts similar in spirit, but differing in the implementation, have been discussed before in the literature basil1 , again resulting in systems of integro-differential equations. In $`\mathrm{\Omega }`$ we define the relations between the potentials $`\varphi `$ and $`\psi `$ and the physical fields as

$$𝐄_\mathrm{\Omega }\varphi _\mathrm{\Omega },𝐃_\mathrm{\Omega }\psi $$
(5)
while in $`\mathrm{\Sigma }`$ we have

$$𝐄_\mathrm{\Sigma }\varphi _{\mathrm{\Sigma }_\mathrm{\Omega }}$$
(6)
while the dielectric displacement field in $`\mathrm{\Sigma }`$ can be represented in terms of a scalar and a vector field arfken ,

$$𝐃_\mathrm{\Sigma }\epsilon _0_\mathrm{\Sigma }𝑑𝐫^{}\epsilon (𝐫,𝐫^{})^{}\varphi _\mathrm{\Sigma }\psi _\mathrm{\Sigma }+\times \xi _\mathrm{\Sigma }.$$
(7)
The scalar field $`\psi `$ thus serves as the potential of the rotation free part $`\stackrel{~}{𝐃}=\psi _\mathrm{\Sigma }`$ of the dielectric displacement field $`𝐃`$. In our setting $`\varphi _\mathrm{\Sigma }`$ is determined by $`\psi _\mathrm{\Sigma }`$ alone, and neither $`\psi _\mathrm{\Omega }`$ nor $`\varphi _\mathrm{\Omega }`$ are affected by $`\xi _\mathrm{\Sigma }`$. Note that $`\xi _\mathrm{\Sigma }`$ can of course be computed from e.g. eq. (7) as soon as $`\varphi _\mathrm{\Sigma }`$ is known. Since we here are interested only in the electrostatic potential $`\varphi `$ and quantities derived from it, we will ignore $`\xi _\mathrm{\Sigma }`$ in the following.

With this definition, the differential equations and boundary conditions, eqs.(2 \- 7) can be brought into the form

$$\mathrm{\Delta }\psi _\mathrm{\Omega }=\varrho ,𝐫ϵ\mathrm{\Omega },\mathrm{\Delta }\psi _\mathrm{\Sigma }=0,𝐫ϵ\mathrm{\Sigma }$$
(8)
and

$$\psi _\mathrm{\Sigma }|_n=\psi _\mathrm{\Omega }|_n,𝐫ϵ\mathrm{\Gamma },\varphi _\mathrm{\Sigma }|_t=\varphi _\mathrm{\Omega }|_t,𝐫ϵ\mathrm{\Gamma }$$
(9)
In addition, there are now two equations relating the potential fields $`\varphi `$ and $`\psi `$,

$$\epsilon _0\epsilon _\mathrm{\Omega }\varphi _\mathrm{\Omega }=\psi _\mathrm{\Omega },𝐫ϵ\mathrm{\Omega }$$
(10)
and eq.(7), to be fulfilled in $`\mathrm{\Sigma }`$.

So far, the introduction of $`\psi `$ has made it possible to rewrite the boundary conditions in a completely local way. The nonlocal problem now only is with eq.(7) in which the nonlocal integral kernel still remains. The latter expression can be simplified further under additional assumptions on the explicit form of $`\epsilon (𝐫,𝐫^{})`$. We now assume that it can be written in the form

$$\epsilon (𝐫,𝐫^{})=\epsilon _{\mathrm{}}\delta (𝐫𝐫^{})+\stackrel{~}{\epsilon }𝒢(𝐫,𝐫^{})$$
(11)
where $`\stackrel{~}{\epsilon }=(\epsilon _\mathrm{\Sigma }\epsilon _{\mathrm{}})/\lambda ^2>0`$, and the Green function $`𝒢`$ solves

$$𝒢(𝐫,𝐫^{})=\delta (𝐫𝐫^{})$$
(12)
for a differential operator $``$ with constant coefficients. Note that $`\epsilon _{\mathrm{}}`$ refers to the value of the dielectric function on smallest scales (i.e., $`𝐫𝐫^{}`$) and can be related to the frequency spectrum of the dielectric function korn4 . While this construction clearly restricts the theory of nonlocal electrostatics to a certain class of dielectric functions, this restriction is not problematic as it applies to most situations of physical interest.

Under these assumptions the application of the differential operator $``$ to eq.(7) yields

$$\epsilon _0\left[\epsilon _{\mathrm{}}\stackrel{~}{\epsilon }\right]\varphi _\mathrm{\Sigma }=\psi _\mathrm{\Sigma }.$$
(13)
As the fields $`\varphi `$ and $`\psi `$ are determined only up to an arbitrary constant, we can drop the gradients on both sides of eq.(13) and are then left with the expression

$$\epsilon _0\left[\epsilon _{\mathrm{}}\stackrel{~}{\epsilon }\right]\varphi _\mathrm{\Sigma }=\psi _\mathrm{\Sigma }.$$
(14)
Eqs.(8 \- 10,14) constitute the novel formulation of nonlocal electrostatics based entirely on partial differential equations.

In order to apply the equations to specific physical situations, the nonlocal dielectric function and thus the Green function of the differential operator $``$ need to be determined. To keep the computations simple – since here we focus only on the basic conceptual features of our approach – we further assume that the dielectric function is isotropic, i.e., $`𝒢(𝐫,𝐫^{})𝒢(𝐫𝐫^{})`$, which is exact far from any boundary. The use of more general expressions is clearly permitted in our theory and is in fact needed for the treatment of realistic situations ritschel ; attard ; korn5 . Within our approach, they lead to more complex Green functions and corresponding differential operators, and will be discussed in a detailed study later.

A standard model for an isotropic nonlocal dielectric function is given by the so-called Fourier-Lorentzian model with a Yukawa-type kernel in real space. The corresponding Green function reads

$$𝒢(𝐫𝐫^{})=\frac{1}{4\pi }\frac{e^{\frac{|𝐫𝐫^{}|}{\lambda }}}{|𝐫𝐫^{}|}$$
(15)
with the differential operator $``$ being given by $`\mathrm{\Delta }1/\lambda ^2`$. With this choice eq.(14) reads

$$\epsilon _0\left[\epsilon _{\mathrm{}}\lambda ^2\mathrm{\Delta }\epsilon _\mathrm{\Sigma }\right]\varphi _\mathrm{\Sigma }=\psi _\mathrm{\Sigma }$$
(16)
Note that due to $`\mathrm{\Delta }\psi _\mathrm{\Sigma }=0`$, no differential operator appears on the rhs of eq.(16).

This result is interesting for two reasons. First, it illustrates that in the limit $`\lambda 0`$, i.e. on length scales large compared to the scale of the orientational correlations, the local limit is recovered. Second, for the differential operator chosen, the form of the equation is apparently that of a Debye-Hückel equation in which the role of the Debye-Hückel screening length is played by the combination $`\lambda (\epsilon _{\mathrm{}}/\epsilon _\mathrm{\Sigma })^{\mathrm{\hspace{0.17em}1}/2}`$, and the potential $`\psi `$ plays the role of the density of mobile charges. We can thus interpret $`\psi `$ as a density of polarization charges in the bulk whose gradient gives rise to the rotation free part of the displacement field $`𝐃`$.

We now turn to illustrate how our formulation of nonlocal electrostatics can be put to use. First, we consider the simplest case of a charge $`q`$ placed at the center of a spherical shell of radius $`a`$. Inside the shell, we assume $`\epsilon _\mathrm{\Omega }=1`$. This system serves as a model for ion solvation korn4 . The equations for $`\psi `$ and $`\varphi `$ can now be solved as follows. In the nonlocal case, the role of the Poisson equation for $`\varphi `$ is taken over by the equation for $`\psi `$ according to eqs.(8,9). The $`\psi `$-potential inside the shell is given by $`\psi _\mathrm{\Omega }=\frac{q}{4\pi r}`$, with the same form given outside, $`\psi _\mathrm{\Sigma }=\frac{q^{}}{4\pi r}`$. The tangential boundary condition is trivially fulfilled, while the normal boundary condition at $`r=a`$ leads to $`q=q^{}`$. Due to the radial symmetry of the problem, eq. (16) reads

$$\epsilon _0\left[\epsilon _{\mathrm{}}\lambda ^2\left(\frac{1}{r^2}\frac{}{_r}\left(r^2\frac{}{_r}\right)\right)\epsilon _\mathrm{\Sigma }\right]\varphi _\mathrm{\Sigma }=\frac{q}{4\pi r}$$
(17)
which is solved by

$$\varphi _\mathrm{\Sigma }(r)=\frac{1}{4\pi \epsilon _\mathrm{\Sigma }\epsilon _0}\frac{q}{r}(1+A\mathrm{exp}(\gamma r))$$
(18)
where the coefficients $`A`$ and $`\gamma `$ follow from coefficient matching and the continuity of $`\varphi `$ at the boundary, $`\varphi _\mathrm{\Sigma }(a)=\varphi _\mathrm{\Omega }(a)`$

$`A{\displaystyle \frac{\epsilon _\mathrm{\Sigma }\epsilon _l}{\epsilon _l}}{\displaystyle \frac{\mathrm{sinh}\nu }{\nu }},\nu \sqrt{{\displaystyle \frac{\epsilon _\mathrm{\Sigma }}{\epsilon _l}}}{\displaystyle \frac{a}{\lambda }},\gamma {\displaystyle \frac{\nu }{a}}=\sqrt{{\displaystyle \frac{\epsilon _\mathrm{\Sigma }}{\epsilon _l}}}{\displaystyle \frac{1}{\lambda }}`$
The electrostatic potential can be used to estimate the solvation energy of monovalent and divalent ions. From $`\varphi _\alpha (r)`$ and $`\psi _\alpha (r)`$, with $`\alpha \{\mathrm{\Omega },\mathrm{\Sigma }\}`$, we can easily compute the free energy of solvation for this setting as the difference of the electrostatic energies in water and vacuum (where the local computation of course remains valid) from $`\mathrm{\Delta }G=\frac{1}{2}\{\rho \varphi \rho \varphi _{\text{vac}}\}𝑑r`$, where the integrals are split into integrals over $`\mathrm{\Omega }`$ and $`\mathrm{\Sigma }`$.

The result of this calculation is shown in Figure 2, where we have compared our results to a corresponding local computation (the Born-model). The correlation length $`\lambda `$ serves as an adjustment parameter; it is the only one in the theory. Like the ion radii, it can also be obtained from microscopic simulations. Our result compares favourably to the experimental data taken from exp . The value for $`\lambda `$ was taken to be $`24.13`$Å; for the ion radii we chose the values according to Åqvist AAqvist . We also tested the set of Shannon radii Shannon1 ; Shannon2  without significant differences on our results. A detailed discussion of the choice of ion radii within nonlocal electrostatics can be found in arxiv .

The simple radially symmetric problem of ion solvation is, of course, not representative for the general character of the solutions of eq.(16), as the tangential boundary condition is trivially fulfilled in this case. In order to elucidate the effect of this boundary condition, we consider the potential generated by a charge $`q`$ placed at a distance $`d`$ from a planar dielectric phase boundary, as sketched in Figure 3.

Within local electrostatics, this problem can be solved by the method of image charges. The symmetry of the situation allows the use of cylindrical coordinates. At a point $`P`$ inside $`\mathrm{\Omega }`$, the local dielectric medium, the electrostatic potential $`\varphi `$ is given by $`\varphi _\mathrm{\Omega }(\rho ,z)=\frac{1}{4\pi \epsilon _0\epsilon _\mathrm{\Omega }}\left(\frac{q}{R_1}+\frac{q^{}}{R_2}\right)`$ where $`R_1=\sqrt{\rho ^2+(dz)^2}`$ and $`R_2=\sqrt{\rho ^2+(d+z)^2}`$ are the lengths of the vectors pointing from the charge $`q`$ located at $`z=d`$, and the image charge $`q^{}`$ at $`z=d`$ to the point $`P`$. Inside the charge-free halfspace $`\mathrm{\Sigma }`$ the potential is given by $`\varphi _\mathrm{\Sigma }=\frac{1}{4\pi \epsilon _0\epsilon _\mathrm{\Sigma }}\frac{q^{\prime \prime }}{R_1}`$ with the image charge $`q^{\prime \prime }`$. The application of the boundary conditions for the normal and tangential components then yields the relations between the image charges $`q+q^{}=\frac{\epsilon _\mathrm{\Omega }}{\epsilon _\mathrm{\Sigma }}q^{\prime \prime },qq^{}=q^{\prime \prime }`$ so that the potential is determined explicitly in both halfspaces. Inside $`\mathrm{\Omega }`$, e.g., one has

$$\varphi _\mathrm{\Omega }(\varrho ,z)=\frac{q}{4\pi \epsilon _0\epsilon _\mathrm{\Omega }}\left(\frac{1}{R_1}+\frac{\epsilon _\mathrm{\Omega }\epsilon _\mathrm{\Sigma }}{\epsilon _\mathrm{\Omega }+\epsilon _\mathrm{\Sigma }}\frac{1}{R_2}\right)$$
(20)
In the nonlocal case, the eqs. for $`\psi `$ have the same form as those for $`\varphi `$ in the local case, so we assume the solutions to be similar. Under this assumption, the normal boundary condition remains unchanged while the tangential boundary condition leads to

$$_\varrho \varphi _\mathrm{\Sigma }=\frac{1}{4\pi \epsilon _0\epsilon _\mathrm{\Omega }}\frac{(q+q^{})}{[\varrho ^2+d^2]^{\mathrm{\hspace{0.17em}3}/2}}\varrho $$
(21)
at $`z=0`$. This equation can be readily integrated along $`\varrho `$ and determines the potential at the dielectric boundary. For $`\varrho ,z\mathrm{}`$, on the other hand, $`\varphi _\mathrm{\Sigma }`$ must vanish. These two conditions then allow to compute the potential $`\varphi _\mathrm{\Sigma }`$ from eq.(16) within $`\mathrm{\Sigma }`$. Since this requires a numerical computation, we leave it for a future publication where we discuss the numerical treatment of our equations. Here, we only give the lowest order effect nonlocality within $`\mathrm{\Sigma }`$ has on $`\mathrm{\Omega }`$. For $`\varrho \sqrt{2}d`$, the electrostatic potential inside $`\mathrm{\Omega }`$ has the same form as in the local theory, but with a ‘renormalized’ bulk value of the dielectric function $`\epsilon _\mathrm{\Sigma }=\epsilon _{\mathrm{\Sigma },loc}2\epsilon _{\mathrm{}}(\lambda ^2/d^2)`$. Transverse variations of the permittivity in the nonlocal medium thus induce a change of the local permittivity proportional to $`(\lambda /d)^2`$ in the vicinity of the boundary, an effect which vanishes for $`\lambda 0`$, and for deeply buried charges, $`d\mathrm{}`$.

As a final remark on the applications of the nonlocal theory of electrostatics we mention the case in which mobile charges are present in the medium surrounding the cavity. Within the linear mean-field theory of local electrostatics, they can be described by Poisson-Boltzmann theory. This stays true within the nonlocal theory presented here. The Boltzmann distribution of the charges simply modifies the rhs of eq.(8). The nonlocal approach can then e.g. be used to quantify recent experimental results of AFM measurements on force-deflection curves at charged mica substrates in water and solutions of monovalent ions Teschke01 . Within a simplified treatment of the dielectric function of water, the orientational effects governed by $`\lambda `$ are found to be on the order of 10 nm Blossey03 . We expect that a more realistic structural model for water, which will become computationally tractable due to our approach, will lead to a much improved description of water orientation near charged surfaces, at least in cases where the assumption of mean-field behaviour is justified.

To conclude, we have presented a novel formulation of nonlocal electrostatics, which includes the effects of spatial dispersion in the dielectric permittivity on surfaces embedded in a solvent, by reformulating it in terms of a two-potential model. While the resulting equations still need to be solved numerically even for simple geometries, this task can now be performed by standard methods developed for partial differential equations. Due to the generality of eq.(14), dielectric functions of greater complexity than the simple radially symmetric choice used here for illustrative purposes can be treated, provided they can be related to known Green functions. Work in this direction is under way.

Acknowledgement. We thank the DFG for support under its research cluster “Informatics methods for the analysis and interpretation of large genomic datasets”, grant LE952/2-3.