format: text

text:
# A High-Frequency and Multi-Epoch VLBI Study of 3C 273

## 1. Introduction

VLBI imaging at millimeter wavelengths (mm-VLBI) gives highest angular resolution and facilitates direct studies of the self-absorbed jet base in AGN. For 3C 273 ($`z=0.158`$) images with an angular resolution of up to $`50`$$`\mu `$as (1 $`\mu `$as $`=10^6`$ arcsec) are made at 86 GHz. This corresponds to a spatial scale of $`1000`$ Schwarzschild radii (for M $`=10^9`$ M$`_{\mathrm{}}`$ ). The small observing beam allows to accurately determine positions of jet components. This facilitates detailed studies of the jet structure and kinematics near the nucleus, in particular with regard to the broad-band (radio to Gamma-ray) flux density variability and the birth of new ‘VLBI components’. Here we summarize new results from our multi epoch (1990 – 1997) and high observing frequency (15, 22, 43 and 86 GHz) VLBI monitoring.

## 2. Results

At sub-milliarcsecond resolution, 3C 273 shows a one sided core-jet structure of several milliarcseconds length. The jet breaks up into multiple VLBI components, which – when represented by Gaussian components – seem to separate at apparent superluminal speeds from the stationary assumed VLBI core. The cross-identification of the model-fit components, seen at different times and epochs, is facilitated by small ($`<0.2`$ mas), and to first order negligible, opacity shifts of the component positions relative to the VLBI-core. Quasi-simultaneous data sets (cf. Fig. 1) demonstrate convincingly the reliability of the component identification, which results in a kinematic scenario, in which all detected jet components (C6 – C18) move steadily (without ‘jumps’ in position) away from the core (Fig. 4). For the components with enough data points at small ($`<2`$ mas) and large ($`>2`$ mas) core separations, quadratic fits to the radial motion r(t) (but also for x(t) –right ascension, and y(t) –declination) represent the observations much better than linear fits. Thus the components seem to accelerate as they move out. The velocities range typically from $`\beta _{app}=48`$ (for $`H_0=100`$ km s-1 Mpc-1, $`q_0=0.5`$).

Dual-frequency maps obtained in 1995 (22/86 GHz) and 1997 (15/86 GHz) allow to measure the spectral index gradient along the jet (see Fig. 2). The spectrum oscillates between $`1.0\alpha +0.5`$ ($`S_\nu \nu ^\alpha `$). Most noteworthy, the spectral gradients did not change significantly over the 2 year time period, although different jet components occupied this jet region. Thus, the geometrical (eg. relativistic aberration) and/or the physical environment (pressure, density, B-field) in the jet must determine the observed properties of the VLBI components. Hence, the latter do not form ‘physical entities’, but seem to react to the physical conditions of the jet fluid.

Further evidence for a fluid dynamical interpretation comes from a study of the mean jet axis and the transverse width of the jet. Both oscillate quasi-sinusoidally on mas-scales. At 15 GHz, the variation of the ridge-line with time could be determined from 4 VLBA maps obtained during 1995 – 1997 (Fig. 3). The maxima and minima of the ridge-line are systematically displaced. This ‘longitudinal’ shift suggests motion with a pattern velocity of $`\beta _{app}=4.2`$. The sinusoidal curvature of the jet axis, however, is more indicative for jet rotation rather than for longitudinal waves. Helical Kelvin-Helmholtz instabilities propagating in the jet sheath could mimic such rotation, which, when seen in projection, would explain also the spectral index oscillation in Figure 2.

For many AGN, a correlation between flux density variability and ejection of VLBI components is suggested. In 3C 273 we identified 13 jet components (C6 – C18) and traced their motion back to their ejection from the VLBI core. The typical measurement uncertainty for the ejection times $`t_0`$ ranges between 0.2 – 0.5 yr. In Figure 5 (left) we plot $`t_0`$ and the millimeter-variability (22 – 230 GHz). We also add the Gamma-ray detections of 3C 273 from EGRET. In Figure 5 (right) we plot the onset times of the mm-flares derived from these light curves (Türler et al. 1999) together with the VLBI ejection time ($`t_0`$) and the Gamma-ray fluxes. For each onset of a mm-flare, we find that a new jet component was ejected.

Although the time sampling of the Gamma-ray data is quite coarse, a relation between component ejection and high Gamma-ray flux appears very likely (note that each Gamma-detection already means higher than usual $`\gamma `$-brightness). From a more detailed analysis (Krichbaum et al. 2001) we obtain for the time lag between component ejection and onset of a mm-flare: $`t_0t_0^{\mathrm{mm}}=0.1\pm 0.2`$ yr. If we assume that the observed peaks in the Gamma-ray light-curve are located near the times $`t_0^\gamma `$ of flux density maxima, we obtain $`t_0^\gamma t_0^{\mathrm{mm}}=0.3\pm 0.3`$ yr. Although the Gamma-ray variability may be faster, this result is fully consistent with the more general finding of enhanced Gamma-ray fluxes mainly during the rising phase of millimeter flares. We therefore suggest the following tentative sequence of events: $`t_0^{\mathrm{mm}}t_0t_0^\gamma `$ – the onset of a millimeter flare is followed by the ejection of a new VLBI component and, either simultaneously or slightly time-delayed, an increase of the Gamma-ray flux. If we focus only on those VLBI components, which were ejected close to the main maxima of the Gamma-ray light-curve in Figure 5, we obtain time lags of $`t_0^\gamma t_0`$ of $`0.5`$ yr for C12, $`0.9`$ yr for C13, $`0.2`$ yr for C16 and $`0.1`$ yr for C18. In all cases the Gamma-rays seem to peak a little later than the time of component ejection. With $`\beta _{app}4`$ near the core, the Gamma-rays would then escape at a radius $`r_\gamma 0.1`$ mas. This corresponds to $`r_\gamma 2000`$ Schwarzschild radii (for a $`10^9`$ M$`_{\mathrm{}}`$ black hole) or $`610^{17}`$ cm, consistent with theoretical expectations, in which Gamma-rays escape the horizon of photon-photon pair production at separations of a few hundred to a few thousand Schwarzschild radii.

References:
Türler M., Courvoisier T.J.-L., & Paltani S., ‘Modelling the submillimeter-to-radio flaring behaviour of 3C 273’, 1999, A&A, 349, p. 45–54.
Krichbaum T.P., et al., ‘High Frequency VLBI Observations of 3C 273’, 2001, A&A, submitted.


text:
# Cosmological Constraints onTerrestrial Planet Formation

## Cosmological Constraints on
Terrestrial Planet Formation

### 0.1 Cosmoplanetology

During a marathon, a runner often does not know which position he is in. His nearest competitors are either so far ahead or so far behind that he doesn’t know where he stands. We earthlings are in a similar situation. We do not yet have a cosmic context for life. We don’t even know if there are other competitors. Approximately 5% of the sun-like stars surveyed possess close-orbiting giant planets (Marcy & Butler 2000) but we still cannot verify if our Solar System is a typical planetary system. Cosmology can help.

Since planets form around stars, the planet formation rate in the universe is strongly linked to the star formation rate. We now have estimates for the star formation rate in the universe (Fig. 1B). If terrestrial planets (‘earths’) always formed around stars then the earth formation rate would be equal to the star formation rate. But earths are made out of metals and metals are the accumulated waste product of stars. Thus the first stars had no earths and the earth formation rate is correlated with the time integral of the star formation rate. The presence of close-orbiting giant planets is both incompatible with the existence of earths and strongly correlated with high metallicity of the host stars (Fig. 2). Thus, there may be a metallicity selection effect: early in the universe with little metallicity, earths are unable to form for lack of material; later on, in star forming regions of very high metallicity, giant planets destroy earths.

The central idea of cosmoplanetology is to piece together a consistent scenario based on current estimates of the star formation rate of the universe, the metallicity evolution of the star-forming regions of the universe and the most recent observations of extrasolar planets. The precision of all of these data sets is improving rapidly, but they can already be combined to yield an estimate of the age distribution of earth-like planets in the universe (Lineweaver 2001). The earth-like planets in the universe are, on average, $`1.8\pm 0.9`$ billion years older than the Earth. If life forms readily on earth-like planets – as suggested by the rapid appearance of life on Earth – this analysis gives us an age distribution for life on such planets and a rare clue about how we compare to other life which may inhabit the universe.

### 0.2 Build it and they will come:
Earth Production $`=`$ Life Production?

The cratering history of the Moon tells us that the Earth underwent an early intense bombardment by planetesimals and comets from its formation 4.56 Gyr ago until $`3.9`$ Gyr ago. For the first 0.5 Gyr, the bombardment was so intense (temperatures so high) that the formation of early life may have been frustrated (Maher & Stevensen 1988). The earliest isotopic evidence for life dates from the end of this heavy bombardment $`3.9`$ billion years ago (Mojzsis et al. 1996). Thus, life on Earth seems to have arisen as soon as temperatures permitted.

Analogous considerations apply to the origin of life in the universe. The oldest isotopic evidence for life in the Universe is $`11`$ Gyr. By ‘oldest isotopic evidence’, I mean, when was the first time in the history of the universe in which the isotopes necessary for the formation of an earth-like planet existed in sufficient abundance to form an earth. The analysis of Lineweaver (2001) indicates that the first earth-like planets probably formed about 2.4 billion years after the big bang in the most metal-rich star forming regions of the universe. It took $`1`$ billion years for the first stars to form and another $`1.4`$ billion years for sufficient metallicity to build up.

To interpret the earth production rate of the universe in Fig. 1B as the life production rate several assumptions need to be made. Among them are: 1) the dominant harbours for life in the universe are on the surfaces of earths in classical habitable zones. 2) life is based on molecular chemistry and cannot be based on just hydrogen and helium. 3) other time-dependent selection effects which promote or hamper the formation of life (supernovae rate?, gamma ray bursts?, cluster environments?) are not as important as the metallicity selection effect discussed here.

In a lottery the odds of winning can be a thousand to one or a million to one. However, no matter what the odds are, the more lottery tickets you buy, the better your chances. Similarly, we do not know how likely life is to form on earth-like planets. It may be very likely or it may be next to impossible. But whatever the odds, the more earths there are, the more likely life will be to form. Thus, the age distribution in Fig. 1B would still represent the age distribtion of life in the universe, independent of how likely such events are.

References

| 1. | Lineweaver, C.H. 1999, Science, 284, 1503 ‘A Younger Age for the Universe’ |
| --- | --- |
| 2. | Lineweaver, C.H. 2001, Icarus in press, astro-ph/0012399 ‘An Estimate of the Age Distribution of Terrestrial Planets in the Universe: Quantifying Metallicity as a Selection Effect’ |
| 3. | Maher, K.A. & Stevenson, D.J. Nature, 331, 612-614 1988 ‘Impact Frustration of the Origin of Life’ |
| 4. | Marcy, G. W. & Butler, R.P. 2000, PASP 112, 137-140 “Millennium Essay: Planets Orbiting Other Suns” |
| 5. | Mojzsis, S.J. et al. 1996, Nature, 384, 55-59 ‘Evidence for life on Earth before 3800 million years ago” |


text:
# Quantum Hall Stripe States in a Tilted Magnetic Field

## 1 Introduction

Recently, a strongly anisotropic transport has been observed in high quality GaAs/AlxGa1-xAs single heterojunctions horst:aps93 ; lilly:prl99a ; du:ssc99 ; pan:prl99 ; lilly:prl99b  at filling factors $`\nu `$ = 9/2, 11/2, etc., and in 2D hole systems shayegan:physicae99  starting at $`\nu =5/2`$. The origin of the magnetotransport anisotropy has not been firmly established yet. The most appealing interpretation suggests that the 2D electron gas spontaneously breaks the translational symmetry by forming a unidirectional charge density wave, as predicted by Hartree-Fock theory koulakov:prl96 ; moessner:prb96 . Because of uncertainty about the reliability of this Hartree-Fock prediction, a special emphasis is placed on tests of its ability to explain experimental results on stripe orientation in tilted magnetic fields.

Many-body RPA/Hartree-Fock calculations using a model, parabolic quantum well system showed jungwirth:prb99  that the stripe orientation relative to the in-plane field is not universal. In narrow quantum wells, stripes are expected to orient perpendicular to the in-plane field. For wider quantum wells, with two occupied subbands at zero tilt angle, a more complex behavior is predicted. Isotropic 2D electron systems with half-filled lowest Landau level (LL) of the second subband can undergo a transition into a UCDW state due to level (anti)crossings induced by the in-plane magnetic field. The orientation of these stripe states can change from parallel to perpendicular as the strength of the in-plane field increases.

In this paper we present calculations of the field-tilt UCDW anisotropy energy that include a self-consistent local-spin-density-approximation (LSDA) description of one-particle states in experimental sample geometries jungwirth:prb99 . Our results are compared to magnetotransport experiments on GaAs/AlxGa1-xAs single heterojunctions with different 2D electron gas densities pan:prl99 ; lilly:prl99b ; lilly:unpubl  and on a square quantum well sample with two occupied electric subbands pan:prl00 .

## 2 Theoretical model and results

The energy per electron of the UCDW state at fractional filling $`\nu ^{}`$ of the valence LL is given by koulakov:prl96

$$E=\frac{1}{2\nu ^{}}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{\Delta }_n^2U\left(\frac{2\pi n}{a}\widehat{e}\right),\mathrm{\Delta }_n=\nu ^{}\frac{\mathrm{sin}(n\nu ^{}\pi )}{n\nu ^{}\pi },$$
(1)
where $`a`$ is the period of the UCDW state, $`\widehat{e}`$ is the direction of charge variation, and $`\mathrm{\Delta }_n`$ is the Fourier transform of the the guiding center occupation function at wave vector $`n2\pi /a`$. In Hartree-Fock theory, the UCDW state energy depends only on $`a`$ and $`\widehat{e}`$ and the optimal UCDW is obtained by minimizing Eq.(1) with respect to these parameters. In Eq.(1), $`U(\stackrel{}{q})`$ can be separated into direct, $`H(\stackrel{}{q})`$, and exchange, $`X(\stackrel{}{q})`$, contributions with

$`H(\stackrel{}{q})`$ $`=`$ $`{\displaystyle \frac{1}{2\pi \mathrm{}^2}}V(\stackrel{}{q})`$
$`X(\stackrel{}{q})`$ $`=`$ $`{\displaystyle \frac{d^2p}{(2\pi )^2}e^{i(p_xq_yp_yq_x)\mathrm{}^2}V(\stackrel{}{p})},`$ (2)
where $`V(\stackrel{}{q})`$ is the RPA-screened effective Coulomb interaction jungwirth:prb99 .

To evaluate $`V(\stackrel{}{q})`$ we choose the in-plane component $`B_{ip}`$ of the magnetic field to be in the $`\widehat{x}`$-direction and use the following Landau gauge for the vector potential, $`\stackrel{}{A}=(0,B_{}xB_{ip}z,0)`$. The one-particle orbitals are calculated using the self-consistent LSDA method and for any sample geometry, with the $`z`$-dependent confining potential, can be written as $`e^{iky}/\sqrt{L_y}\phi _{i,\sigma }(x\mathrm{}^2k,z)`$, where $`k`$ is the wave vector which labels states within LL $`i`$, $`\sigma `$ is the spin index, and $`\mathrm{}^2=\mathrm{}c/eB_{}`$. The translational symmetry responsible for LL degeneracy leads to a 2D wavefunction $`\phi _{i,\sigma }(x,z)`$ which is independent of the state label $`k`$, except for the rigid shift by $`\mathrm{}^2k`$ along $`x`$-axis. This in turn leads to two-particle matrix elements of the Coulomb interactions with a dependence on state labels which is identical to that for the lowest LL of a zero-thickness 2D electron system provided the 2D Coulomb interaction is replaced by the effective interaction $`V(\stackrel{}{q})`$.

For an infinitely narrow electron layer the orbital effects of $`B_{ip}`$ vanish and the 2D Coulomb interaction, $`V(\stackrel{}{q})`$, reduces to $`e^{q^2\mathrm{}^2/2}/q`$ $`(L_N(q^2\mathrm{}^2/2))^2`$ $`2\pi e^2\mathrm{}/ϵ`$ where $`L_N(x)`$ is the Laguerre polynomial and $`N`$ is the orbit radius quantum number. Starting from $`N=1`$, zeros of $`L_N(q^2\mathrm{}^2/2)`$ occur at finite $`q=q^{}`$, producing a zero in the repulsive Hartree interaction at wave vectors where the attractive exchange interaction is strong. For the half-filled valence LL the corresponding UCDW state consists of alternating occupied and empty stripes of electron guiding center states with a modulation period $`a2\pi /q^{}`$.

In finite-thickness 2D systems the dependence of the effective interaction on wavevector magnitude $`q`$ and orientation $`\varphi `$ relative to the in-plane field direction can be accurately approximated jungwirth:prb99  by $`V(\stackrel{}{q})=V_0(q)+V_2(q)`$ $`\mathrm{cos}(2\varphi )`$. At weak $`B_{ip}`$, the isotropic term $`V_0(q)`$ has a wavevector-dependence similar to that of the effective interaction in the infinitely narrow 2D layer. The corresponding curve for the valence LL at $`\nu =9/2`$ for a single heterojunction with 2D density $`n_e=2.67\times 10^{11}`$ cm-2 lilly:prl99b  is shown in the lower inset of Figure 1.

The non-zero anisotropy coefficient $`V_2(q)`$ of the effective interaction at $`B_{ip}>0`$, also shown in the inset, is responsible for the non-zero UCDW anisotropy energy $`E_A`$, defined jungwirth:prb99  as the total Hartree-Fock energy of stripes oriented parallel with $`B_{ip}`$ minus the total energy of stripes perpendicular to $`B_{ip}`$. The main plot of Figure 1 shows $`E_A`$ calculated for single heterojunction samples. The two curves correspond to high and low density 2D electron systems with one occupied electric subband in perpendicular magnetic field. The subband wavefunctions are plotted in the upper inset. We found that stripes are oriented perpendicular to the in-plane field for the considered interval of densities $`n_e=1.52.67\times 10^{11}`$ cm-2, consistent with available experimental data on single heterojunctions pan:prl99 ; lilly:prl99b ; lilly:unpubl . Similar agreement between theory and experiment is also abtained for other half-filled high LLs.

In Figure 2 we present results for a 350 $`\AA `$ wide GaAs quantum well with $`n_e=4.6\times 10^{11}`$ cm-2 pan:prl00 . At $`\nu =9/2`$ and 11/2, and $`B_{ip}=0`$, the valence LL is the $`N=0`$ LL of the second subband. The effective interaction $`V_0(q)`$ has no zeros at finite wavevector, as seen in the inset of Figure 2, and the UCDW state is not expected to form, consistent with experiment pan:prl00 . An in-plane field $`B_{ip}3`$ T induces a level anticrossing accompanied by a development of a minimum in $`V_0(q)`$ at finite $`q`$. As discussed above, it is the minimum of the interaction energy at finite wavevector that opens the possibility for the formation of the UCDW state. The theoretical and experimental pan:prl00  in-plane fields corresponding to the onset of the UCDW and transport anisotropy, respectively, are remarkably close.

The main graph in Figure 2 indicates that, in this sample, the sign of $`E_A`$ changes at $`B_{ip}10`$ T. The stripes orient parallel with $`B_{ip}`$ at low fields while perpendicular orientation is obtained in the high in-plane field region, consistent with the experimentally observed interchange of easy and hard current axes at similar $`B_{ip}`$.

We conclude that the close agreement between complex experimental data and theoretical results for all studied samples leave little doubt as to the UCDW origin of the observed transport anisotropies.


text:
# Experimental evidence of chaotic advection in a convective flow

## I Figure captions


text:
# The Los Alamos Spallation Driven Solid Deuterium Ultracold Neutron Source: Working Notes of Mar. 30, 1998 with Embellishment

## I Prologue

The idea that solid deuterium could be a useful source of ultracold neutron (UCN) dates to the early 1980’s when R. Golub et al.  proposed the “thin film” UCN source. The basic idea is that a material, such as solid deuterium (SD2) has a finite neutron absorption cross section, but a high rate of UCN production compared to superfluid helium. A small amount of SD2 contained in a UCN storage cell and irradiated by cold neutrons will produce UCN until the rate of production equals the rate of absorption. Because both processes are proportional to the amount of SD2, the steady-state UCN density is independent of the amount of material. Although the production rate is higher for SD2 than for superfluid helium, the ultimate density of UCN in a superfluid helium filled long lifetime storage cell is around 10-100 times greater.

However, there are applications where SD2 will be useful for a UCN converter. It has been suggested that the thin film source would be easier to use near the core of a high-flux reactor, primarily because less material is required to achieve a particular UCN current.

A specific application of SD2 to UCN production at a spallation neutron source was brought to our attention in Jan. 1998 by A. Serebrov  at a workshop held near St. Petersburg, Russia. The basic idea is that the UCN storage bottle need only be connected to the SD2 during a short time after the spallation pulse; one can then enjoy a high production rate together with an effectively increased storage lifetime. Upon returning to Los Alamos we (D. Bowman, S. Lamoreaux, C. Morris) discussed the “UCN Factory” described in  and came to the conclusion that it had a number of problems that include wasting phase space of the cold neutron flux, the necessity of a UCN window between the SD2 and storage vessel, the use of an unwieldy amount of SD2 (approaching 1 m3), an overly complicated shutter mechanism, and it seemed that the SD2 thickness was such that UCN might not be able to propagate out of the material irrespective of the window. Such a system might eventually be made to work, but without understanding the basic physical mechanisms involved, it seemed unwise to invest the large amount of resources required to build even a scaled version of the Factory.

In addressing these problems, in particular the cold neutron phase space, we were challenged as to why they were indeed problems. In answering the challenge, the following model source was proposed. Upon further reflection, the geometry is close to ideal for our specific application.

Our model source closely follows a proposal by Pokotilovski  who really came up with the idea of separating the production and storage regions, with a shutter that opens only around the time of the spallation pulse. R. Golub had provided S. Lamoreaux with a copy of this paper in 1995, but we had no conscious memory of this paper at the time of the St. Petersburg workshop; the paper emerged when we were digging through our notes. Our model source has only minor embellishments over the original Pokotilovski proposal which the interested reader can find by comparing these notes to .

## II Notes of Mar. 20, 1998

### A Basic scheme of source

A rough schematic of the physics demonstration UCN source is shown in Fig. 1.

### B Cold neutron production:

15 neutrons/800 MeV proton from spallation,

$$P_c=15\frac{I}{e}10^{14}/\mathrm{sec}\mathrm{at}1\mu \mathrm{A}$$
(1)
($`P_c`$ is the cold neutron production rate).

If the poly or LD2 moderator is 30 K or lower, the produced neutrons will moderate to $``$30 K (moderation processes turn off for neutron temperatures below $``$25-30 K).

For poly, the albedo is about 90% for 30 K neutrons. We can estimate the mean free path of a cold neutron in the flux trap region as $``$ 40 cm before loss (up the UCN guide) or absorption. This gives a cold neutron dwell time in the flux trap region

$$\tau _c=\frac{0.4\mathrm{m}}{700\mathrm{m}/\mathrm{s}}=0.66\mathrm{ms}.$$
(2)
The cold neutron density in the flux trap region is them

$$\frac{P_c\tau _c}{V_c}=\rho _c.$$
(3)
Take $`V_c`$ as 1.5 times solid D2 volume $`V_D`$

$$V_c=1.5\mathrm{}.$$
(4)
So, at 1 $`\mu `$A,

$$\frac{10^{14}}{\mathrm{sec}}0.66\times 10^3\mathrm{sec}\times \frac{1}{1500\mathrm{c}\mathrm{c}}=4.4\times 10^7/\mathrm{cc}$$
(5)
Then, the cold neutron flux is

$$\frac{1}{4}\rho _cv_c;\mathrm{with}\mathrm{v}_\mathrm{c}=700\mathrm{m}/\mathrm{s}.$$
(6)
$$\mathrm{\Phi }_0=8\times 10^{11}\mathrm{n}/\mathrm{cm}^2\mathrm{s}$$
(7)
### C UCN Production rate

See Fig 3.17 of Ultracold Neutrons (Golub, Richardson, Lamoreaux) . For the thin film source, $`\rho =150\times 10^{11}\mathrm{\Phi }_0/\mathrm{cm}^3`$.

The production rate (of UCN) is $`P_u=\rho /\tau `$ and $`\tau =0.12\mathrm{sec}`$ at $`T<4`$ K. So, $`P_u=1.25\times 10^8\mathrm{\Phi }_0`$ for SD2 in a Be bottle. We must correct for the potential; $`U_{SD_2}=108`$ neV, this is compensated by the vertical rise. For the thin film source, $`P(U_{Be}^{3/2}U_{SD_2}^{3/2})`$ for our source, $`P^{}U_{Fomblin}^{3/2}`$

$$\frac{P}{P^{}}=\frac{252^{3/2}108^{3/2}}{160^{3/2}}=1.42$$
(8)
so four our system $`P_u=9\times 10^9\mathrm{\Phi }_0`$.

Total number of UCN produced is

$$V_DP_u=10009\times 10^98\times 10^{11}=7.2\times 10^6\mathrm{UCN}/\mathrm{sec}.$$
(9)
### D UCN lifetime and UCN density

We have

$$\frac{7.2\times 10^6\mathrm{UCN}/\mathrm{sec}}{15000\mathrm{cc}\mathrm{total}}=480\frac{\mathrm{UCN}}{\mathrm{sec}\mathrm{cm}^3}.$$
(10)
The lifetime is given approximately by the ratio of the total volume $`V_T`$ to the solid D2 volume, times the UCN lifetime in solid D2.

$$\frac{V_T}{V_D}\tau _D=\frac{15}{1}0.12\mathrm{sec}=1.8\mathrm{sec}.$$
(11)
This gives

$$\rho _{ucn}=\frac{480\mathrm{UCN}}{\mathrm{cm}^3\mathrm{s}}1.8\mathrm{s}=764\mathrm{UCN}/\mathrm{cm}^3.$$
(12)
### E Equilibration time etc.

The equilibrium time for the system is 1.8 sec, the shortest lifetime in the system. Thus, the 1 $`\mu `$A beam would have to be turned on for 3.6 sec (two lifetimes) after which the shutter to the storage volume is closed.

After 3.6 sec of filling, the storage volume will be filled with $`10^7`$ UCN. If these are released at a constant rate over 400 sec, loss will be minimized, and the UCN current will be $`10^7/400=2.5\times 10^4`$/sec.

The 800 MeV beam would have to be applied every 400 sec; this rate could be increaed for a higher UCN flux. Also, the 800 MeV beam current could be increased; it’s a question of radiation damage and heat load. (Power for 1 $`\mu `$A, 800 MeV=800 Watts.)

## III Epilogue

The prototype source was constructed and tested in Aug. 1998. As anticipated, the window created a significant UCN loss; e.g., there is an effectively small shunt impedance in parallel with the source, as shown in Fig. 2.

A windowless system was designed and built; yet, a substantially lower UCN density compared to these notes was obtained. This was due to a number of causes.

The state of affairs as of late 1998 are given in .

First, the flux trap is rather more complicated that the simple estimate above would indicate. The cold flux density extends further into the moderator and Be reflector, so the mean free path is modified. Monte Carlo calculations by R. Hill indicate that the flux trap lifetime is closer to 160 $`\mu `$ sec, and the effective volume is larger. This results in a factor of 10 reduction in cold neutron flux.

In a test of the production rate of UCN by scattering cold neutrons in solid deuterium in March 1998, done at the Hahn-Meitner Institut in Berlin, the rate was within a factor of 3-5 of the expected rate given in Sec. II.C above, assuming all detected neutrons were UCN. However, the expected temperature dependence to the apparent production rate was not observed (possibly due to loss of UCN before they can propagate out of the solid). The observed UCN flux and lack of temperature dependence could be explained by assuming a UCN lifetime of a few milliseconds in the solid. The same lack of temperature dependence was initially observed in the demonstration source; we eventually found that the para molecular state, which persists from its room temperature concentration through freezing at 4 K, can upscatter UCN with a lifetime of about 4 milliseconds .

A para/ortho converter was installed on the demonstration source and the observed UCN production rate agrees well with Monte Carlo calculations, and a temperature dependence is observed.

The hydrogen contamination of available deuterium gas (0.1 at%), upscattering by phonons at 5 K (lowest practical temperature), and 2% para contamination (lowest practical conversion temperature is 17 K at which temperature there is still a reasonable vapor pressure), gives a net lifetime of 25 milliseconds.

A final weakness in the working notes is the assumption that one liter of solid D2 could be used; in fact, the mean free path of a UCN is about 8 cm due to the incoherent cross section alone (and is possibly shorter for an amorphous material) and limits the volume; roughly 200 cc is a practical upper limit for the solid volume for the geometry shown in Fig. 1.

These factors taken together give a factor of 10-100 less UCN than the simple calculation. Nonetheless, the concept still has potential to be the basis of the world’s most intense UCN source, and with minor modifications we expect 300 UCN/cc densities in bottles of volumes 10’s of liters; the practical limit to a spallation pulse for a target assembly amenable to this project is 50 $`\mu `$C in one second; this would be repeated every 10 sec so the average current would be about 5 $`\mu `$A. Questions relating to heating of the sample and generation of the ortho state due to radiation remain to be studied. Details of the Los Alamos effort will soon be published by the collaboration.


text:
# Single hole spectral function and spin-charge separation in the 𝑡-𝐽 model

## Abstract

Worm algorithm Monte Carlo simulations of the hole Green function with subsequent spectral analysis were performed for $`0.1J/t0.4`$ on lattices with up to $`L\times L=32\times 32`$ sites at temperature as low as $`T=J/40`$, and present, apparently, the hole spectral function in the thermodynamic limit. Spectral analysis reveals a $`\delta `$-function-sharp quasiparticle peak at the lower edge of the spectrum which is incompatible with the power-law singularity and thus rules out the possibility of spin-charge separation in this parameter range. Spectral continuum features two peaks separated by a gap $`4÷5t`$.

For almost four decades the problem of hole dynamics in magnetic systems has attracted constant interest with applications ranging from properties of charge carriers in magnetic semiconductors and insulators  to vacancies in solid 3He . The research in this area exploded with the discovery of high temperature superconductors in cuprates, where superconductivity appears upon light doping of AFM insulators. Despite an enormous theoretical effort over the years and quite a variety of treatments (for reviews, see, e.g., Ref. ) a complete solution of this inherently strong-coupling problem still does not exist, especially in the most interesting region of $`t>J`$, where $`J`$ is the exchange coupling constant and $`t`$ is the hopping matrix element in the $`tJ`$ Hamiltonian

$$H=t\underset{<ij>,s}{}c_{is}^{}c_{js}+J\underset{<ij>}{}\left(𝐬_i𝐬_j\frac{1}{4}n_in_j\right).$$
(1)
Here $`c_{j\sigma }`$ is projected (to avoid double occupancy) fermion annihilation operator, $`n_i=_sc_{is}^{}c_{is}2`$ is the occupation number, $`𝐬_i=_{ss^{}}c_{is}^{}\sigma _{ss^{}}c_{is^{}}`$ is spin-$`1/2`$ operator, and $`<ij>`$ denote nearest neighbor sites of the 2D square lattice.

The central problem in the hole dynamics is whether or not its spin and charge degrees of freedom separate. The standard way to answer this question is to study the spectral function $`A_𝐩(\omega )=\pi ^1ImG_𝐩(\omega )`$, where $`G_𝐩(\omega )`$ is the hole Green function. If there is an elementary excitation associated with the hole, the spectral function is supposed to feature a peak at the lower edge of the spectrum. What is crucial, however, is not the presence of the peak itself, but its functional form . Within the self-consistent Born approximation scheme (SCBA)  one finds finite overlap between the bare hole and low-energy quasiparticle states, which means that the peak is $`\delta `$-functional and the hole is described as coherently propagating spin-polaron in the nearly ordered antiferromagnetic (AFM) background (with vanishing scattering at low temperature due to small density of spin waves) .

In contrast to that, various resonating-valence-bond (RVB) descriptions and Anderson’s general arguments about breakdown of the Fermi-liquid picture in the system with no-double-occupancy constraint (see, e.g., Ref. ) strongly suggest that power-law singularity, which is indicative of spin-change separation scenario, might be the case (there is even a claim that the quasiparticle weight $`Z`$ should be rigorously zero ). To make the issue more confusing, angle resolved photoemission spectroscopy experiments  show very broad maximum in $`A(\omega )`$ which can be considered both as the quasiparticle peak with anomalously large broadening or as the evidence for composite nature of quasiparticles .

The quasiparticle picture was supported by exact calculations on small clusters , but system sizes (up to $`32`$ sites) were too small to perform finite-size scaling. Variational calculations, Green function Monte Carlo and density matrix renormalization group studies were mostly concerned with the dispersion law $`ϵ_𝐤`$ (lowest energy in a given momentum sector). Large scale simulations of the imaginary time Green function $`G_𝐤(\tau )`$ were performed recently using a combination of the loop-cluster Monte Carlo method for the AFM state and hole evolution in the fixed space-time spin background . This method works only for relatively large exchange $`J>0.6t`$, since magnetic background is simulated without the hole and polaron-type distortions have to be accounted for as quantum fluctuations before the hole is introduced. For $`J/t<0.6`$ the error bars in $`G_𝐤(\tau )`$ are too large for reliable spectral analysis (see below).

To summarize, we still lack evidence that for small $`J`$ the quasiparticle weight remains finite in the thermodynamic limit and the lowest peak has nothing to do with the power-law singularity. We thus find it important to rule out the possibility that t-J model may explain the data of Refs.  (as suggested by Ref. ), so that extensions of the model such as $`t^{}`$ and $`t^{\prime \prime }`$ terms  or frustrating exchange couplings are proven necessary.

Speaking classically, moving hole breaks AFM bonds and thus its energy increases linearly with the travel distance  (this consideration, or the string-potential picture, is most appropriate for the $`tJ_z`$ model). It is believed that the ground energy scaling $`E_{𝐤_0}J^{2/3}`$ \[where $`𝐤_0=(\pi /2,\pi /2)`$\] and excitation spectrum are described by the string-potential picture , and transverse spin fluctuations do not “erase” strings completely. In the limit of small $`J`$ the theory predicts that several resonances in $`A(\omega )`$ have to be seen with the peak positions being strictly related to the eigenvalue properties of Airy functions. Early exact diagonalization studies on clusters $`4\times 4`$  attributed two peaks above the ground state to string resonances, however later studies on larger clusters  were not able to detect the second resonance, and the spectral function was showing strong size dependence. What happens at small $`J`$ in the thermodynamic limit remains an open question.

In this letter we present results for $`G_𝐤(\tau )`$ and $`A_𝐤(\omega )`$ obtained from Monte Carlo simulations on systems with $`16\times 16`$ and $`32\times 32`$ sites and at temperatures as low as $`T/J=0.025`$ (for the largest system size) using continuous-time Worm algorithm  in combination with the recently developed spectral analysis which is capable of resolving infinitely sharp features in $`A(\omega )`$ . The method itself is free from any systematic errors, and we were unable to detect finite-size corrections in our data; thus, we believe, our results describe correctly the thermodynamic limit. In the parameter range studied $`0.1J/t0.4`$, the lowest peak in $`A_{𝐤_0}(\omega )`$ is a $`\delta `$-function within the resolution limit of order of $`0.01÷0.03t`$, which means that our quasiparticle is the spin-polaron. For the excitation spectrum we observe two well-separated peaks for all values of $`J`$. The ground state energy scaling does follow the $`J^{2/3}`$ law predicted by the string-potential theory. Although we were unable to resolve individual string resonances, we believe that their combined effect is seen as the first peak in the spectral continuum since its position also scales as $`J^{2/3}`$. The high-energy peak is roughly at a constant distance $`5t`$ from the ground state.

Worm algorithm is based on the idea that world-line configurations of spins and the hole are updated through the space-time motion of the creation and annihilation operators. In Fig. 1 and Fig. 2 we show the typical configuration of the Heisenberg AFM with the hole, and Monte Carlo updates which apply to it. Physical configurations contributing to the hole Green function are those which have no spin end points (denoted by filled circles).

Since the world-line representation is based on the expansion of the statistical evolution operator $`e^{H/T}`$ in powers of $`t`$ and $`J`$ it suffers from the sign problem which first appears in order $`t^2J^3`$ (see Fig. 3). It is worth noting that if not for the sign problem, spin-charge separation can be ruled out by the analysis of world-line configurations. Let the hole be created by $`c_i`$ operator. If spin charge separation does take place, one should see, following the evolution of the system configuration in imaginary time, an extra spin density leaving the hole creation site and going into the bulk. i.e., the world-line density far from the hole should increase. At $`T=0`$ the AFM ground state is ordered (as opposite to the spin liquid state) and the minimal possible change in the world-line density is equivalent to having exactly one extra world

line, which can be traced out and interpreted as spin-$`1`$ magnon excitation. We may now construct an operator which has finite overlap with the quasiparticle excitation as a product $`s^{}c_{}`$, where $`s^{}`$ is added to cancel the extra magnon in the bulk. However, up to a single hopping transition the above composite operator is identical to $`c_{}`$, and we conclude that holes are good quasiparticles in contradiction with original assumption (probably rephrasing the proof of Ref. ).

However, in the presence of the sign problem the above consideration should be taken with extreme caution. It may turn out that ordered world-line configurations compensate each other completely and single-configuration conclusions are misleading, as suggested in Ref.  where hole-related sign problem is called the “irreparable phase string effect” and argued to cause spin-charge separation.

The sign problem implies that we may not calculate $`G(\tau )`$ reliably over long time scales and have to restrict our simulation to $`t\tau <3÷4`$ to suppress sign fluctuations by larger statistics. Fortunately, on this time scale $`G_𝐤(\tau )`$ is already in its asymptotic regime and the data are sufficiently accurate to reveal the ground state properties. Formally, calculations are done at finite $`T`$ but its value is more than an order of magnitude smaller than the energy of the lowest magnon state in a given system size. For each value of J the calculation time was about 2 weeks on a PIII-600 workstation.

In Fig. 4 we show simulated $`G_{𝐤_0}(\tau )`$ for $`J/t=0.4`$ and the asymptotic law $`Z_{𝐤_0}e^{E_{𝐤_0}\tau }`$ with the quasiparticle weight and ground state energy obtained from the weight and position of the $`\delta `$-peak in $`A_{𝐤_0}(\omega )`$. Note, that for small values of $`J`$ the data have to be very accurate to describe correctly how $`G(\tau )`$ approaches its asymptotic behavior $`G_{𝐤_0}Z_{𝐤_0}e^{E_{𝐤_0}\tau }`$. Error bars are shown but are smaller than the symbol size (the relative accuracy is better than $`10^2`$ even for points with the largest $`\tau `$ where the sign problem was the most severe).

The spectral analysis of $`G_{𝐤_0}(\tau )`$ was done using stochastic optimization procedure developed earlier for the polaron problem. $`A(\omega )=N^1_{i=1}^NA_i(\omega )`$ is obtained as an average over spectral densities which optimize deviations between $`G(\tau )`$ and $`𝑑\omega e^{\omega \tau }A_i(\omega )`$. The parameter space of $`A_i(\omega )`$ is defined by the step-wise constant functions, which, in particular, includes infinitely sharp peaks and is not associated with any predefined set of frequencies  (it is known that $`\delta `$-peaks can not be handled satisfactorily by the maximum entropy method ).

In Fig. 5 we show our results for $`A_{𝐤_0}(\omega )`$ calculated at points $`J/t=0.4,0.2,0.1`$. We clearly see a $`\delta `$-sharp peak at the lower edge of the spectrum. The structure of this peak is incompatible with the power law singularity since its width is smaller than the lowest magnon excitation in our system \[for $`J/t=0.4`$ the quasiparticle peak width is only $`0.01t`$ (!) while the natural scale for the power law is set by $`J`$\]. This is the central result of our paper which conclusively rules out spin-charge separation scenario for the single hole dynamics in the $`tJ`$ model and confirms finite quasiparticle weight in the thermodynamic limit. To verify that finite-size and finite-temperature corrections are negligible we performed long time simulations in a $`32\times 32`$ lattice at temperature $`T=J/40`$ for $`J/t=0.2`$, but within the error bars $`G(\tau )`$ was indistinguishable from the result obtained for $`L=16`$ and $`T=J/20`$.

Unfortunately, the ill-defined problem of numeric analytic continuation does not allow to study fine structures in the spectral density, especially if they are “screened” by low- and high-frequency peaks. (The low-frequency

peak is fixed by the asymptotic long-time behavior of $`G(\tau )`$, while the high-frequency peak is fixed by the short-time decay of $`G(\tau )`$ where the data are extremely accurate.) Our tests show that multiple peaks in the middle can not be resolved by spectral analysis even when we use analytically exact $`G(\tau )`$ data. It means that the absence of multiple string resonances above the ground state in our results for $`A(\omega )`$ may not  be considered as a proof that string potential picture fails in quantum case. We would rather consider the second peak as a “course grain” description of spectral density at intermediate energies. However, if string excitations do exist, their combined effect should be seen as the $`J^{2/3}`$ scaling law for the peak position. In Fig. 6 we plot peak positions as functions of $`(J/t)^{2/3}`$ for $`0.1J/t1.2`$ with error bars obtained as peak half-widths. We conclude that for the second peak the scaling law is obeyed within the error bars. The high-energy peak stays roughly at a constant distance from the ground state, and clearly the physics behind it is different.

We thank O. Ruebenacker, P. Stamp, and V. Kashurnikov for valuable discussions. This work was supported by the National Science Foundation under Grant DMR-0071767.


text:
# Reply to Comment

In the comment  Eggeling, Vollbrecht and Wolf suspect our method in Ref.  is not practical. Here we explain our result and method and show that our example can tell one how to judge a separable state, and so our method is practical, at least for many mixed states.

A. Explanation

Recently we considered the necessary and sufficient condition of separability and gained a practical criterion . Our main idea is as followings:

1. We decompose a mixed state $`\rho `$ of $`mn`$ dimension Hilbert space into (m-1)(n-1) pairs which are in 2$`2`$ dimension Hilbert space, from Ref., we get the necessary condition of each pair being parallel, i.e. $`\rho `$ being separable.

2. Each parallel pair corresponding to a set of matrices U, if there are intersection of all matrices U, $`\rho `$ is separable, otherwise $`\rho `$ is inseparable. We shown by our example that judging the intersection existing or not comes down to judging a set of special symmetric quadratic equations having solutions or not, which is not difficult, so our method is practical.

B. Reply

The first step above generalizes Wootter’s work to higher dimension. This is meaningful not only for separable criterion but also for other approach(we are working in this). Parallel pair is different from PPT, though they are equivalent for 2$`2`$ system, but in the case of higher dimension they have different character and usage.

If one has found a matrix U which brings up a set of separable pure states decomposition of $`\rho ,`$one can say $`\rho `$ is separable. But how to judge this matrix existing or not? The second step above has told one how to find this matrix U, so it is practical.

We think if the authors of the comment consider the example in Ref carefully, they will find our method is practical.


text:
# Synthesis and Processing of 𝑀𝑔𝐵₂ powders and wires

## I Introduction

The discovery of superconductivity in the intermetallic compound, $`MgB_2`$ near $`40K`$ has prompted both fundamental studies of the mechanism causing superconductivity and practical studies of supercurrent transport at grain boundaries. Bud’ko and co-workers have shown the existence of a boron isotope effect in that $`T_cM^\alpha `$ with $`\alpha =0.26`$, where $`T_c`$ is the superconducting transition temperature and $`M`$ is the isotopic mass. This result implies that the electron-phonon interaction is important in the determination of $`T_c`$ and that light elements with high phonon frequencies might be important in the search for new metallic superconductors with high transition temperatures. Energy gap measurements and neutron scattering measurements of the phonon spectrum all point to a conventional mechanism for the superconductivity. On the practical side, Canfield and co-workers have shown a method to convert commercially available boron fibers into very low resistivity $`MgB_2`$ wire and Takano and co-workers have measured the current carrying capacity of hot-pressed pellets of the material. Other measurements of flux pinning and intergrain coupling have been made using studies of flux creep.

In this paper, we report a systematic study of the vapor diffusion method of sample preparation for both microcrystalline powders and boron ($`B`$) fibers, which are available commercially in kilometer lengths. The primary variable in the study is time, and for the wire samples the effect of defects in the starting $`B`$ fibers. Of particular interest is the role of the tungsten-boride core in creating cracks and transmitting $`Mg`$ along the length of the fiber.

## II Experiment

For powder synthesis, we place a stoichiometric mixture of bright lumps of $`Mg`$ ($`99.9\%`$) and fine crystalline $`{}_{}{}^{11}B`$ powder ($`99.5\%`$ ) that is less than $`100\mu m`$ in diameter into a $`Ta`$ tube whose ends are sealed under a partial atmosphere of $`Ar`$. For the wire synthesis, we cut $`B`$ fiber lengths several $`cm`$ long and either $`100\mu m`$ or $`141\mu m`$ in diameter and place them in $`Ta`$ tubes under a partial pressure of $`Ar`$ with excess $`Mg`$. These starting $`B`$ fibers have a tungsten-boride core about $`15\mu m`$ in diameter with breaks every few $`cm`$.

After the $`Ta`$ tubes are sealed, they are in turn sealed in quartz ampoules with $`175mbar`$ of $`Ar`$ and placed in a box furnace preheated to $`950^{}C`$. It is believed that using a preheated furnace inhibits the growth of higher boron phases. After reaction, the samples are rapidly cooled by putting the quartz tube in running cold water. The formation of $`MgB_2`$ from $`B`$ causes the powders to expand and bow the $`Ta`$ tube outward and causes $`100\mu m`$ diameter wires to expand by more than $`30\%`$ in diameter.

Samples are characterized by X-ray diffraction, magnetization, and electrical resistivity. X-ray diffraction spectra are taken with $`CuK_\alpha `$ radiation in a Scintag diffractometer where the scan angle $`2\theta `$ varies from $`20^{}`$ to $`90^{}`$with a scan step of $`0.020^{}`$. Magnetization is measured in a Quantum Design $`MPMS5`$ magnetometer with a $`6cm`$ sample travel in a magnetic field of $`25Oe`$ in a zero-field-cooled mode. Resistivity is measured with a standard four probe technique using an $`LR700`$ resistance bridge.

## III Results and discussion

The lattice constants of the $`MgB_2`$ phase, determined from X-ray data, change relatively little with time of reaction for powders and have values close to those previously reported. The dynamics of the growth of the $`MgB_2`$ phase are illustrated by the study of powder diffraction X-ray data after different times of reaction in Fig. 1. The major $`B`$ peak at $`38^{}`$, marked by the solid line on the figure, is small and goes away during processing. For particle sizes used here (less than $`100\mu m`$), the $`MgB_2`$ peaks grow quickly between $`30min`$ and $`60min`$ and remain relatively unchanged from $`120min`$ to $`240min`$. The $`Si`$ peaks used for a standard have been removed from these spectra. It should be noted that when the reaction ampoule was placed in a furnace, heated to $`950^{}C`$, left for $`120min`$, and quenched, the powder X-ray spectra showed many second phase peaks. Zero-field-cooled magnetization measurements for these same samples show the transition width narrowing from $`15min`$ to $`30min`$ and then remaining relatively constant at about $`0.5K`$.

Electrical resistivity of the sintered pellets is compared to the wires in Fig.2. The inset in Fig. 2 shows an expanded view near the transition temperature and indicates values at $`40K`$ of about $`1.5\mu ohmcm`$ for the pellets and $`0.4\mu ohmcm`$ for the wire segments. These processing methods give very low electrical resistivity in the normal state for both pellets and wires. Because both the relatively open structure of the pellet and the dense structure of the wire samples have comparable residual resistivity ratios, $`RRR20`$, we think the lower resistivity values for the wire sample are due to its higher density.

In order to study the synthesis of $`MgB_2`$ wire, a series of exposure time, t, of the boron fibers in $`Mg`$ vapor was measured for $`15mint240min`$. For the processing of the $`100\mu m`$ diameter fibers, the scanning electron micrographs in Fig. 3a and 3b illustrate the behavior seen. In these micrographs, $`Mg`$ has a much higher atomic number than $`B`$ so the shading is a good indication of $`Mg`$ content. In Fig. 3a, a length of fiber reacted for $`60min`$ has been polished through a section that misses the tungsten-boride core in the middle of the fiber. The $`Mg`$ concentration is very high at the outside surface region marked (1), and there is a clear front of $`Mg`$ moving across the fiber. At the boundary between the regions marked (2) and (3), energy dispersive X-ray spectra (EDS) show that the $`Mg`$ level abruptly falls by about $`40\%`$. Region (3) has a fairly high $`Mg`$ content, even though it has not transformed to the $`MgB_2`$ phase. As shown by the diameter lines, the diameter of this section has grown by more than $`10\%`$ from the original $`100\mu m`$.

A different section of fiber reacted for $`60min`$ in the same run as in Fig. 3a is shown in Fig. 3b. Here the polishing just barely exposes the tungsten-boride core along the centerline. This SEM micrograph shows many of the common defects. The dark line along the center is a void region, which often occurs near the tungsten-boride core, and an EDS spectrum shows very strong $`W`$ lines in this region. The bright region marked (A) is an $`MgB_2`$ region and the dark region marked (B) is an unreacted $`B`$ region. If optical pictures are taken, the $`B`$ regions stand up above the $`MgB_2`$ regions because they are harder and polish less rapidly. Figure 3c shows an EDS line scan measurement of the $`Mg`$ concentration taken along the dark line running from upper left to lower right in the left-center portion of Fig. 3b. The $`Mg`$ count on Fig. 3c drops from about $`1.4`$ to about $`0.9`$ going from the bright region to the dark region. The tick mark on the line in Fig. 3b corresponds to the vertical line in Fig. 3c and to a point where the scan moves from a dark region back into a bright region near the core. Magnesium seems to diffuse very rapidly along the core region and diffuse outward from there to form $`MgB_2`$. The void region shown by the dark area in the upper left corner of Fig. 3b is unusually large, and the bottom of the void is clearly seen in the SEM. Voids are commonly about $`1020\mu m`$ across, and the bottom of the void generally can be seen in the SEM. The diagonal regions of high $`Mg`$ content we believe are cracks where the $`Mg`$ can diffuse rapidly and start the formation of $`MgB_2`$.

A $`B`$ fiber reacted for $`30min`$ is shown in Fig. 4a. Again the polishing depth does not cut through the core. The sample has much less conversion to $`MgB_2`$, shown by the bright spots. EDS measurements show the highest $`Mg`$ content in the bright spot marked (C), substantially less $`Mg`$ near the edge at the point marked (A), and even less at the point near the center marked (B).

Another $`B`$ fiber reacted for $`15min`$ is shown in Fig. 4b. In the left edge of this fiber, the tungsten core shows through, and the diameter of the fiber is still close to $`100\mu m`$. The EDS data show only a small amount of $`Mg`$ near the edge of the fiber and practically none in the center of the fiber.

When exposure time, $`t=120min`$, the $`100\mu m`$ starting diameter fibers are fully reacted and end up having a final diameter of about $`160\mu m`$. To examine how the process scales with initial fiber size, a $`141\mu m`$ diameter fiber was exposed for a much longer time, $`t=36h`$, with the results shown in Fig. 5. The conversion to $`MgB_2`$ has increased the diameter of the wire to greater than $`180\mu m`$, as shown by the diameter measurements. There are no second phase particles large enough to be seen in the SEM micrograph, and the fiber is relatively void free.

## IV Conclusions

Low resistivity, high purity powders and highly dense wires of $`MgB_2`$ can be synthesized by exposure of boron powder and boron fibers to $`Mg`$ vapor at $`950^{}C`$ for adequate exposure time. Comparable results are seen for both $`100\mu m`$ and $`141\mu m`$ diameter fibers. For $`B`$ fibers that are placed in a furnace preheated to $`950^{}C`$, the $`MgB_2`$ phase begins to form quickly after about $`15min`$ and is about $`50`$ percent complete after an hour. Diffusion of $`Mg`$ along the tungsten-boride core and along voids or cracks in the fiber are important means of transport. As the $`MgB_2`$ phase forms, the diameter of the fiber grows from about $`141\mu m`$ to about $`190\mu m`$. As reported earlier, the wire segments provide a highly conducting material with a normal state resistivity of about $`0.4\mu ohmcm`$ and a superconducting critical current density of about $`20,000A/cm^2`$ at $`1T`$ and $`20K`$.

## V Acknowledgments

Work is supported by the U.S. Department of Energy, Basic Energy Sciences, Office of Science, through the Ames Laboratory under Contract No. W-7405-Eng-82.


text:
# The missing stress-geometry equation in granular media

## 1 Introduction

Granular media are ubiquitous yet complex materials with puzzling properties . The simplest model of a static granular material is that where grains are considered to be perfectly hard, perefectly rough and each grain $`\alpha `$ has a coordination number $`z^\alpha =d+1`$ (where $`d`$ is the dimension of the system). Under these conditions Newton’s equations of intergranular force and couple balance can be solved . The system is in the state of mechanical equilibrium and particles can not experience deformation under load so there is no displacement field present. Thus the only immediately available macroscopic equation has the form

$$_j\sigma _{ij}(𝐫)=g_i(𝐫),$$
(1)
where $`\sigma _{ij}(𝐫)`$ is the macroscopic stress tensor and $`g_i(𝐫)`$ is external force at the boundaries. The vector equation (1) gives $`d`$ equations for $`\frac{d(d+1)}{2}`$ components of $`\sigma _{ij}(𝐫)`$ leaving $`\frac{d(d1)}{2}`$ further equations required to solve for the macroscopic stress tensor. Thus one should be able to derive them from the geometry of the contact network which is assumed to be specified. In order to put the above comments into formulae we draw a diagram (see Figure 1) of one grain in contact with 4 nearest neighbours (i.e. grains that are in contact with the reference grain $`\alpha `$).

The geometrical specification of the system is given by the set of $`\frac{_\alpha ^Nz^\alpha d}{2}`$ contact points $`\{𝐑^{\alpha \beta }\}`$. The centroid of contacts of the reference grain $`\alpha `$ is defined by vector $`𝐑^\alpha `$

$$𝐑^\alpha =\frac{_\beta 𝐑^{\alpha \beta }}{z^\alpha },$$
(2)
where the summation sign $`_\beta `$ means the sum over all nearest neighbours of the grain $`\alpha `$. The distance between grains $`\alpha `$ and $`\beta `$ is defined as the distance between their centroids of contacts

$$𝐑^{\alpha \beta }=𝐑^\beta 𝐑^\alpha =𝐫^{\alpha \beta }𝐫^{\beta \alpha },$$
(3)
where $`𝐫^{\alpha \beta }`$ is the vector joining the centroid of contact with the contact point. The second vector which characterises the relative position of neighbouring centroid with respect to the contact point is defined by (see Figure 2)

$$𝐐^{\alpha \beta }=(𝐫^{\alpha \beta }+𝐫^{\beta \alpha }).$$
(4)
Newton’s laws of force and couple balance for every grain give us the system of $`\frac{Nd(d+1)}{2}`$ equations for $`\frac{zdN}{2}`$ interparticle forces $`𝐟^{\alpha \beta }`$ (see Figure 3)

$$\underset{\beta }{}f_i^{\alpha \beta }=g_i^\alpha ,$$
(5)
$$f_i^{\alpha \beta }+f_i^{\beta \alpha }=\mathrm{\hspace{0.17em}0},$$
(6)
$$\underset{\beta }{}ϵ_{ikl}f_k^{\alpha \beta }r_l^{\alpha \beta }=c_i^\alpha .$$
(7)
where $`g_i^\alpha `$ is the external force acting on grain $`\alpha `$ and $`c_i^\alpha `$ is the external couple.

The tensorial force moment $`S_{ij}^\alpha `$ for grain $`\alpha `$ is defined as

$$S_{ij}^\alpha =\underset{\beta }{}f_i^{\alpha \beta }r_j^{\alpha \beta },$$
(8)
and it is symmetric tensor given that $`c_i^\alpha =0`$. The macroscopic stress tensor $`\sigma _{ij}(𝐫)`$ can be obtained by averaging $`S_{ij}^\alpha `$ over the packing

$$\sigma _{ij}(𝐫)=S_{ij}^\alpha ,$$
(9)
where $`\mathrm{}=\frac{1}{V}_{\alpha =1}^N\mathrm{}\delta (𝐫𝐑^\alpha )`$ in the simplest case and $`V`$ is the volume of the packing. The method offered by the authors in  was to consider the probability functional for the set $`\{S_{ij}^\alpha \}`$

$$\begin{array}{cc}\hfill P\left\{S_{ij}^\alpha \right\}& =\underset{\alpha ,\beta }{}\delta \left(S_{ij}^\alpha \underset{\beta }{}f_i^{\alpha \beta }r_j^{\alpha \beta }\right)P\left\{f^{\alpha \beta }\right\}𝒟f^{\alpha \beta },\hfill \end{array}$$
(10)
where

$$\begin{array}{cc}\hfill P\left\{f^{\alpha \beta }\right\}& =𝒩\underset{\alpha ,\beta }{}\delta \left(\underset{\beta }{}f_i^{\alpha \beta }g_i^\alpha \right)\hfill \\ & \times \delta \left(\underset{\beta }{}ϵ_{ikl}f_k^{\alpha \beta }r_l^{\alpha \beta }\right)\hfill \\ & \times \delta \left(f_i^{\alpha \beta }+f_i^{\beta \alpha }\right),\hfill \end{array}$$
(11)
and the normalisations, $`𝒩`$ and $``$ are functions of contact network configuration. The main goal of Ref. was to transform (10) into

$$\begin{array}{c}\hfill P\left\{S_{ij}^\alpha \right\}=\underset{\alpha }{}\delta \left(\underset{\beta }{}K_{ijk}^{\alpha \beta }S_{jk}^\beta g_i^\alpha \right)\delta \left(\underset{\beta }{}P_{ijkl}^{\alpha \beta }S_{kl}^\beta \right)\end{array}$$
(12)
where delta-functions contain the complete system of equations for the set of tensorial force moments $`\{S_{ij}^\alpha \}`$. The method of Ref. was to exponentiate all delta-functions in (10)

$$\begin{array}{c}\hfill P\left\{S_{ij}^\alpha \right\}=\underset{\alpha ,\beta }{\overset{N}{}}e^{iA}𝒟f^{\alpha \beta }𝒟\zeta ^\alpha 𝒟\gamma ^\alpha 𝒟\lambda ^\alpha 𝒟\eta ^{\alpha \beta },\end{array}$$
(13)
where

$$\begin{array}{cc}\hfill A& =\underset{\alpha }{}\zeta _{ij}^\alpha \left(S_{ij}^\alpha \underset{\beta }{}f_i^{\alpha \beta }r_j^{\alpha \beta }\right)\hfill \\ & +\gamma _i^\alpha \left(\underset{\beta }{}f_i^{\alpha \beta }g_i^\alpha \right)\hfill \\ & +\lambda _i^\alpha \left(\underset{\beta }{}ϵ_{ikl}f_k^{\alpha \beta }r_l^{\alpha \beta }\right)\hfill \\ & +\eta _i^{\alpha \beta }\left(f_i^{\alpha \beta }+f_i^{\beta \alpha }\right).\hfill \end{array}$$
(14)
After integrating out the fields $`f_i^{\alpha \beta }`$, $`\lambda _i^\alpha `$ and $`\eta _i^{\alpha \beta }`$ the following system of linear equations for the conjugate fields $`\{\zeta _{ij}^\alpha \}`$ and $`\{\gamma _i^\alpha \}`$ was obtained

$$\zeta _{ij}^\alpha r_j^{\alpha \beta }\gamma _i^\alpha =\zeta _{ij}^\beta r_j^{\beta \alpha }\gamma _i^\beta .$$
(15)
From (15) it was shown in Ref. that $`\zeta _{ij}^\alpha `$ has the representation

$$\zeta _{ij}^\alpha =\zeta _{ij}^{\alpha \mathrm{\hspace{0.17em}0}}+\zeta _{ij}^\alpha ,$$
(16)
where a particular solution $`\zeta _{ij}^{\alpha \mathrm{\hspace{0.17em}0}}`$ gave the first delta-function in (12)

$$\underset{\beta }{}S_{ij}^\alpha M_{jl}^\alpha R_l^{\alpha \beta }\underset{\beta }{}S_{ij}^\beta M_{jl}^\beta R_l^{\beta \alpha }=g_i^\alpha ,$$
(17)
and a complimentary function $`\zeta _{ij}^\alpha `$ which satisfies the following system of linear equations

$$\zeta _{ij}^\alpha r_j^{\alpha \beta }\zeta _{ij}^\beta r_j^{\beta \alpha }=0.$$
(18)
This system of linear equations for $`\{\zeta _{ij}^\alpha \}`$ gave the required $`\frac{Nd(d1)}{2}`$ constraints on $`\{S_{ij}^\alpha \}`$. The present paper concentrates on $`\zeta _{ij}^\alpha `$ and uses the result of Ref. that

$$P\left\{S_{ij}^\alpha \right\}=\underset{\alpha =1}{\overset{N}{}}\delta \left(\underset{\beta }{}S_{ij}^\alpha M_{jl}^\alpha R_l^{\alpha \beta }S_{ij}^\beta M_{jl}^\beta R_l^{\beta \alpha }g_i^\alpha \right)P\left\{S_{ij}^\alpha \right|geometry\},$$
(19)
where $`P\left\{S_{ij}^\alpha \right|geometry\}`$ contains the set of $`\frac{d(d1)}{2}`$ missing equations which in continuum limit might take a form

$$P_{ijkl}\sigma _{kl}+_jT_{ijkl}\sigma _{kl}+_j_lU_{ijkl}\sigma _{km}+\mathrm{}=0.$$
(20)
In the first term we have $`P_{ijkl}=P_{jikl}`$ which gives the correct number of equations. The missing equation (20) is of purely geometric origin. We deliberately avoid using the term ”constitutive relation” (for there is no deformation or displacement in this model) and call Eq. (20) the stress-geometry equation. The authors believe this kind of situation offeres a new kind of challenge in theoretical physics and although this paper presents a solution, it does not have the elegance and completeness that one might expect from the solution derived by means of some variational principle. There may be some analogy here to the problems of dynamics where, although the Lagrange equations with Lagrange multipliers will solve any non-holonomic problem, the really powerful way is to use the Gibbs-Appell equations .

## 2 The missing stress-geometry equation

In order to obtain $`P\left\{S_{ij}^\alpha \right|geometry\}`$ and derive (20) we need to solve Eq.(18). Let us describe the method of Ref.. We notice that in order to obtain the precise number of missing equations (which is $`\frac{Nd(d1)}{2}`$), Eq.(18) appears to be too many equations. Because $`\{𝐫^{\alpha \beta }\}`$ satisfy the linear relation $`_\beta 𝐫^{\alpha \beta }=0`$ from Eq.(2) there are many internal identities and careful counting shows that it can only contain $`Nd`$ equations. For example if when solved Eq. (18) were to give

$`\zeta _{11}^\alpha +\zeta _{22}^\alpha =\mathrm{\hspace{0.17em}0}`$ (21)
$`\zeta _{12}^\alpha =\mathrm{\hspace{0.17em}0}`$
then

$$P\left\{S_{ij}^\alpha \right|geometry\}=\underset{\alpha }{\overset{N}{}}e^{i_\alpha ^N(S_{11}^\alpha \zeta _{11}^\alpha +\mathrm{\hspace{0.17em}2}S_{12}^\alpha \zeta _{12}^\alpha +S_{22}^\alpha \zeta _{22}^\alpha )}𝒟\zeta ^\alpha ,$$
(22)
gives $`S_{11}^\alpha =S_{22}^\alpha `$ and no constraint on $`S_{12}^\alpha `$. One can force Eq.(18) into only two equations (for $`d=2`$) by taking scalar product with vectors $`a_i^{\alpha \beta }`$ and $`b_i^{\alpha \beta }`$ and then summing over $`\beta `$:

$$\underset{\beta }{}a_i^{\alpha \beta }\zeta _{ij}^\alpha r_j^{\alpha \beta }\underset{\beta }{}a_i^{\alpha \beta }\zeta _{ij}^\beta r_j^{\beta \alpha }=\mathrm{\hspace{0.17em}0},$$
(23)
$$\underset{\beta }{}b_i^{\alpha \beta }\zeta _{ij}^\alpha r_j^{\alpha \beta }\underset{\beta }{}b_i^{\alpha \beta }\zeta _{ij}^\beta r_j^{\beta \alpha }=\mathrm{\hspace{0.17em}0},$$
(24)
which gives us $`2N`$ equations. Suppose that we regard the differences between $`\zeta _{ij}^\alpha `$ and $`\zeta _{ij}^\beta `$ to be expandable in $`𝐑^{\alpha \beta }`$.

$$\zeta _{ij}^\beta =\zeta _{ij}^\alpha +R_k^{\alpha \beta }\frac{\zeta _{ij}^\alpha }{R_k^\alpha }+\mathrm{}.$$
(25)
The first approximation is then

$$\underset{\beta }{}a_i^{\alpha \beta }\zeta _{ij}^\alpha R_j^{\alpha \beta }=\mathrm{\hspace{0.17em}0},$$
(26)
and similarly from Eq.(24)

$$\underset{\beta }{}b_i^{\alpha \beta }\zeta _{ij}^\alpha R_j^{\alpha \beta }=\mathrm{\hspace{0.17em}0}.$$
(27)
We have declared in Ref. that there are two obvious vectors to be used in Eqs.(26) and (27), namely $`𝐑^{\alpha \beta }`$ and $`𝐐^{\alpha \beta }`$. Hence we obtain two configuration tensors (analogous, but different from the fabric tensors used in the soil mechanics literature )

$$F_{ij}^\alpha =\underset{\beta }{}R_i^{\alpha \beta }R_j^{\alpha \beta },$$
(28)
and

$$G_{ij}^\alpha =\frac{1}{2}\left(\underset{\beta }{}Q_i^{\alpha \beta }R_j^{\alpha \beta }+Q_j^{\alpha \beta }R_i^{\alpha \beta }\right),$$
(29)
then we have

$$\zeta _{ij}^\alpha F_{ij}^\alpha =\mathrm{\hspace{0.17em}0},$$
(30)
$$\zeta _{ij}^\alpha G_{ij}^\alpha =\mathrm{\hspace{0.17em}0},$$
(31)
and after exponentiating Eqs.(30,31) in $`P\left\{S_{ij}^\alpha \right|geometry\}`$ and eliminating the auxilirary fields one finds the missing stress-geometry equation

$$\left|\begin{array}{ccc}S_{11}^\alpha & F_{11}^\alpha & G_{11}^\alpha \\ \\ S_{12}^\alpha & F_{12}^\alpha & G_{12}^\alpha \\ \\ S_{22}^\alpha & F_{22}^\alpha & G_{22}^\alpha \end{array}\right|=0.$$
(32)
The first approximation (26,27) (christened the ”first coordination shell approximation”) can be illustrated in the following way: we rearrange Eq.(18) and after multiplying by $`𝐑^{\alpha \beta }`$ and $`𝐐^{\alpha \beta }`$ obtain

$$\zeta _{ij}^\alpha \underset{\beta }{}R_i^{\alpha \beta }R_j^{\alpha \beta }+\underset{\beta }{}(\zeta _{ij}^\alpha \zeta _{ij}^\beta )r_j^{\beta \alpha }R_i^{\alpha \beta }=0,$$
(33)
$$\zeta _{ij}^\alpha \underset{\beta }{}Q_i^{\alpha \beta }R_j^{\alpha \beta }+\underset{\beta }{}(\zeta _{ij}^\alpha \zeta _{ij}^\beta )r_j^{\beta \alpha }Q_i^{\alpha \beta }=0,$$
(34)
It is clear that if the second term in Eqs.(33) and (34) is neglected one obtains Eqs.(30) and (31). A naive attempt to transform (32) into the macroscopic equation for $`\sigma _{ij}(𝐫)`$ by averaging $`S_{ij}^\alpha `$, $`F_{ij}^\alpha `$ and $`G_{ij}^\alpha `$ fails (in the case of an isotropic configuration) because

$$G_{ij}^\alpha =\frac{1}{V}\underset{\alpha =1}{\overset{N}{}}\underset{\beta }{}Q_i^{\alpha \beta }R_j^{\alpha \beta }=\frac{1}{V}\underset{\beta =1}{\overset{N}{}}\underset{\alpha }{}Q_i^{\beta \alpha }R_j^{\beta \alpha }=\mathrm{\hspace{0.17em}0}$$
(35)
Thus for an isotropic packing the first term in Eq. (20) vanishes in the first coordination shell approximation. This gives rise to the conditional probability distribution functions. Thus if we are given $`S_{12}^\alpha `$, the probability of finding $`S_{11}^\alpha S_{22}^\alpha `$ is

$$P\left\{S_{11}^\alpha S_{22}^\alpha \right|S_{12}^\alpha \}=\frac{2}{\pi }\frac{|S_{12}^\alpha |}{(S_{11}^\alpha S_{22}^\alpha )^2+(S_{12}^\alpha )^2}.$$
(36)
and vice versa. This distribution can then be introduced for corresponding components of the macroscopic stress tensor subject to the absence of mesoscopic correlations in the packing. For an anisotropic configuration (32) characterized by the distribution of $`\{G_{ij}^\alpha \}`$ with some nonvanishing average yields macroscopic equation when

$$\sigma _{11}\sigma _{22}=\mathrm{\hspace{0.17em}2}\sigma _{12}\text{tan}\varphi .$$
(37)
where $`\varphi `$ is the angle of repose in the case when configuration is prepared in the form of a sandpile . The obvious criticism of the derivation method of Ref. is that one could employ some different vector to obtain Eqs.(30) and (31). For example instead of $`𝐑^{\alpha \beta }`$ one could use $`𝐑^{\alpha \beta }\kappa (𝐑^{\alpha \beta },𝐐^{\alpha \beta })`$ where $`\kappa `$ is any scalar function of $`𝐑^{\alpha \beta }`$ and $`𝐐^{\alpha \beta }`$; or one could go to the next coordination shell of the reference grain $`\alpha `$ and employ $`𝐑^{\beta \gamma }`$ and $`𝐐^{\beta \gamma }`$ (where $`\gamma `$’s are the other two neighbours of $`\beta `$, see Figure 4) and so on. Thus we have offered a path to the missing equation which works also in 3-D, but it is not unique . Presumably the internal symmetries of (18) will lead to the same macroscopic equation when Eqs.(26,27) are used for any $`a_i^{\alpha \beta }`$, $`b_i^{\alpha \beta }`$, but it is not easy to see how. In the next section we will offer some new viewpoints which approach the problem from a different standpoint, but which will confirm the earlier results, and suggest new approaches.

## 3 New methods of derivation

In the previous section we have shown that it is possible to go from Newton’s equations (5-7) to equations for conjugate fields (15) that can be then used to derive the complete set of equations for $`\{S_{ij}^\alpha \}`$. We can reverse the process and from

$$P\left\{S_{ij}^\alpha \right|geometry\}=\underset{\alpha ,\beta }{\overset{N}{}}e^{i_\alpha ^NS_{ij}^\alpha \zeta _{ij}^\alpha }\delta (\zeta _{ij}^\alpha r_j^{\alpha \beta }\zeta _{ij}^\beta r_j^{\beta \alpha })𝒟\zeta ^\alpha ,$$
(38)
obtain

$$P\left\{S_{ij}^\alpha \right|geometry\}=\underset{\alpha ,\beta }{}\delta (S_{ij}^\alpha \underset{\beta }{}p_i^{\alpha \beta }r_j^{\alpha \beta })\delta (p_i^{\alpha \beta }+p_i^{\beta \alpha })𝒟𝐩^{\alpha \beta }.$$
(39)
Since the only constraint on $`\left\{𝐩^{\alpha \beta }\right\}`$ is given by

$$p_i^{\alpha \beta }+p_i^{\beta \alpha }=\mathrm{\hspace{0.17em}0}$$
(40)
the missing stress-geometry equation comes from the condition on $`S_{ij}^\alpha `$ that

$$S_{ij}^\alpha =\frac{1}{2}\underset{\beta }{}\left(p_i^{\alpha \beta }r_j^{\alpha \beta }+p_j^{\alpha \beta }r_i^{\alpha \beta }\right)$$
(41)
can be satisfied by a set of ”pseudo-forces” $`\{𝐩^{\alpha \beta }\}`$ which obey (40), but are free from the constraint of Newton’s second law (5). ”Pseudo-force” $`𝐩^{\alpha \beta }`$ can be written in the following form

$$p_i^{\alpha \beta }=(\psi ^\alpha +\psi ^\beta )R_i^{\alpha \beta }+(\chi ^\alpha \chi ^\beta )Q_i^{\alpha \beta },$$
(42)
which satisfies (40) and provides Eq. (32). We can think of $`\{S_{ij}^\alpha \}`$ as a $`\frac{Nd(d+1)}{2}`$ component vector $`𝐒`$, and $`\{𝐩^{\alpha \beta }\}`$ as a $`\frac{Nd(d+1)}{2}`$ component vector $`𝐏`$. Then the set of equations (41) can be written as

$$𝐒=𝐑𝐏$$
(43)
where $`𝐑`$ is a $`\frac{Nd(d+1)}{2}\times \frac{Nd(d+1)}{2}`$ matrix with a large repetition of elements since $`_\beta 𝐫^{\alpha \beta }=0`$. This matrix has $`\frac{Nd(d1)}{2}`$ zero eigenvalues which gives the right number of constraints on $`\{S_{ij}^\alpha \}`$. The argument here is of the ”must be so” type. Apart from the easy proof that $`\text{Det}𝐑=0`$ (by adding rows) we have not succeeded in proving that the $`\frac{Nd(d+1)}{2}\times \frac{Nd(d+1)}{2}`$ matrix has $`\frac{Nd(d1)}{2}`$ zero eigenvalues, but it ”must be so”! Another and perhaps easier method is to return to the Eq.(18) for $`\zeta _{ij}^\alpha `$ and sum it over $`\beta `$ using $`_\beta 𝐫^{\alpha \beta }=0`$. Therefore we have

$$\underset{\beta }{}\zeta _{ij}^\beta r_j^{\beta \alpha }=\mathrm{\hspace{0.17em}0}.$$
(44)
This can be used in Eq.(39)

$$P\left\{S_{ij}^\alpha \right|geometry\}=\underset{\alpha ,\beta }{\overset{N}{}}e^{i_\alpha ^NS_{ij}^\alpha \zeta _{ij}^\alpha }\delta (\underset{\beta }{}\zeta _{ij}^\beta r_j^{\beta \alpha })𝒟\zeta ^\alpha ,$$
(45)
as in the derivation of (41), but now losing a positional index

$$S_{ij}^\alpha =\frac{1}{2}\underset{\beta }{}\left(\varphi _i^\beta r_j^{\alpha \beta }+\varphi _j^\beta r_i^{\alpha \beta }\right).$$
(46)
As before this equation implies a relationship between the components of $`S_{ij}^\alpha `$. Suppose now that we are looking for an ansatz for $`\varphi _i^\beta `$. Since $`\varphi _i^\beta `$ is a vector it can be represented as a superposition of two obvious candidates $`𝐑^{\alpha \beta }`$ and $`𝐐^{\alpha \beta }`$

$$\varphi _i^\beta =\psi ^\alpha R_i^{\alpha \beta }+\chi ^\alpha Q_i^{\alpha \beta }$$
(47)
where new quantities $`\psi ^\alpha `$ and $`\chi ^\alpha `$ are scalars. This gives us

$$S_{ij}^\alpha =\psi ^\alpha \left(F_{ij}^\alpha +G_{ij}^\alpha \right)+\chi ^\alpha \left(H_{ij}^\alpha +G_{ij}^\alpha \right),$$
(48)
where $`H_{ij}^\alpha =_\beta Q_i^{\alpha \beta }Q_j^{\alpha \beta }`$. Eliminating $`\psi ^\alpha `$ and $`\chi ^\alpha `$ leads to

$$\left|\begin{array}{ccc}S_{11}^\alpha & F_{11}^\alpha +G_{11}^\alpha & H_{11}^\alpha +G_{11}^\alpha \\ \\ S_{12}^\alpha & F_{12}^\alpha +G_{12}^\alpha & H_{12}^\alpha +G_{12}^\alpha \\ \\ S_{22}^\alpha & F_{22}^\alpha +G_{22}^\alpha & H_{22}^\alpha +G_{22}^\alpha \end{array}\right|=0.$$
(49)
Note that $`H_{ij}^\alpha `$ will on average be a multiple of $`F_{ij}^\alpha `$ and this gives us Eq. (32). It is possible to obtain from (46)the third term in (20). The second term is non-vanishing only in special cases of periodic arrays . Let us construct the following interpolation of $`\varphi _i^\beta `$

$$\varphi _i^\beta =\varphi _i^\alpha +R_j^{\alpha \beta }_j\varphi _i^\alpha ,$$
(50)
after substituiting it into (46) and summing it over $`\beta `$ using $`_\beta 𝐫^{\alpha \beta }=0`$ we have

$$S_{ij}^\alpha =_k\varphi _i^\alpha \underset{\beta }{}R_k^{\alpha \beta }r_j^{\alpha \beta }.$$
(51)
Crude averaging gives us

$$\sigma _{ij}=\frac{1}{V}\underset{\alpha }{\overset{N}{}}\underset{\beta }{}r_j^{\alpha \beta }R_k^{\alpha \beta }_k\varphi _i^\alpha =F_{jk}_k\varphi _i.$$
(52)
We can now eliminate $`\varphi _i`$ and obtain the well-known Navier equation which imposes kinematic compatibility on the stresses

$$\frac{^2\sigma _{xx}}{y^2}+\frac{^2\sigma _{yy}}{x^2}\mathrm{\hspace{0.17em}2}\frac{^2\sigma _{xy}}{xy}=\mathrm{\hspace{0.17em}0}.$$
(53)
This equation corresponds to the third term in (20) in the case of an isotropic packing and implies that the stress tensor components can be expressed in terms of the Airy function  whose discrete analogues are $`\psi ^\alpha `$ and $`\chi ^\alpha `$ in (47).This also means that the set of microscopic constraints (46) is consistent with the macroscopic equation (1).

However in 3-D we have encountered mathematical difficulties with this simple averaging procedure. This highlights the puzzle we face. One could put much more complicated versions of $`\varphi _i^\beta `$, for example any functional form which employs vectors constructed out of the vectors of the Figure 2 or indeed the next coordination shell (see Figure 4). The problem is a geometric puzzle. A set of contact points $`\{𝐑^{\alpha \beta }\}`$ corresponds to the packing of grains with coordination number $`z=d+1`$. From these points associated with the reference grain $`\alpha `$ one can construct the vectors of Figure 2, and find the basic equations using $`\zeta _{ij}^\alpha `$ or $`𝐩^{\alpha \beta }`$, or $`\varphi _i^\alpha `$. With these equations one may take ad hoc steps to obtain the complete set of equations for the macroscopic stress tensor $`\sigma _{ij}(𝐫)`$, but we have failed to find a systematic procedure such as is possible for the corrections of the stress-force equation as outlined in Ref.. The problem is now purely geometric, but we emphasize that although real granular materials have many features omitted here, we are studying the simplest possible case and the geometric puzzle offered here although difficult is quite basic, and no simpler case can be found.

## 4 Acknowledgments

We wish to acknowledge the financial support of Leverhulme Foundation S.F.E), ROPA grant from EPSRC (UK) and Research Fellowship from Wolfson College (D. V. G.). Authors thank Prof. R. C. Ball and Dr. R. Blumenfeld (who have derived a different route to solve this problem ) for stimulating discussions.