text:
# Experimental feasibility of measuring the gravitational redshift of light using dispersion in optical fibers

## I INTRODUCTION

Einstein’s Equivalence Principle (EEP) is the pillar underlying all metric theories of gravity, including the theory of general relativity. This principle states that effects of a uniform gravitational field are identical to the effects of uniform acceleration of a coordinate system. EEP can be decomposed into three components: the weak equivalence principle, local Lorentz invariance, and local position invariance (LPI). Precision experiments testing each of these components have yielded results consistent with EEP and the theory of general relativity . However, it is interesting to note that current solar and atmospheric neutrino data are consistent with a small violation of EEP . Given the fundamental nature of EEP, higher precision tests are desirable. Such experiments might rule out some metric theories of gravity, uncover a deviation from our expectation that would require a modification of our current view of gravity, and place constraints on new macroscopic forces (which might be perverse enough to respect the weak equivalence principle, but not LPI).

Local position invariance requires the results of local nongravitational experiments to be independent of position and time. Limits on violations of LPI have been set by searching for variation in the ‘constants’ of nature and by gravitational redshift experiments . To date, all experimental results are consistent with LPI.

In a uniform gravitational field with strength g, the frequency shift of light traveling a distance z along the direction of maximum field gradient is given by

$$\nu ^{}\nu \left\{1\pm \frac{gz}{c^2}\right\}$$
(1)
to first order, where the negative sign is used for light traveling against the gravitational field and the positive sign is used for light traveling with the field. Experiments that measure the gravitational frequency shift of light report results as a limit on the deviation observed from the equation above. The deviation, $`\alpha `$, is defined by

$$\frac{\mathrm{\Delta }\nu }{\nu }=(1+\alpha )\frac{\mathrm{\Delta }U}{c^2},$$
(2)
where $`\mathrm{\Delta }`$U is the change in gravitational potential (gz in the example above) encountered by the light. Currently, the best limit on $`\alpha `$ comes from a 1976 experiment that compared the frequency of a precision hydrogen maser frequency standard on the surface of the earth with one launched by a Scout Rocket to a height of approximately 10,000 km. This experiment yielded a result of $`\alpha `$\<2x10-4 .

## II A new class of experiments

Continuous wave light will travel a vertical distance h in an optical fiber with index of refraction n in a time

$$t=\frac{h}{c}n$$
(3)
The frequency dependence of this flight time is given by

$$dt=\frac{h}{c}\left(\frac{dn}{d\nu }\right)d\nu .$$
(4)
Assuming this to take place near the earth’s surface (let g be constant), the change in the gravitational potential over a distance dh induces a frequency change

$$d\nu =\frac{g\nu }{c^2}dh.$$
(5)
Over a height h, the flight time of the light will differ from that given by equation (3) by an amount

$$\delta t=\frac{g\nu h^2}{2c^3}\left(\frac{dn}{d\nu }\right),$$
(6)
Assuming a sinusoidal form for the light waves traveling down the fiber optic cable of

$$e^{i(\beta z\omega t)},$$
(7)
where

$$\beta =\frac{n\omega }{c}.$$
(8)
and

$$\omega =2\pi \nu .$$
(9)
Defining the phase as

$$\varphi =\beta (\omega )z\omega t,$$
(10)
the phase variation becomes

$$\delta \varphi =(\delta \beta )z+\beta \delta z(\delta \omega )t\omega \delta t.$$
(11)
Specializing to our case, z become the height difference, h, and since this is a constant in our proposal, the variation of h is zero and the second term in the above equation disappears. Now the propagation parameter, $`\beta `$ can be expanded, to first order as

$$\beta (\omega )=\beta _0(\omega _0)+\beta _1(\omega _0)\mathrm{\Delta }\omega +\mathrm{}..$$
(12)
Keeping the first order term in the variation and substituting this and our previously derived equations for the t’s into the phase equation above, we get

$$\delta \varphi =\beta _1\delta \omega h\frac{nh}{c}\delta \omega \omega \frac{g\omega h^2}{2c^3}\frac{dn}{d\omega }$$
(13)
From the definitions above

$$\beta _1=\frac{d\beta }{d\omega }=\frac{\omega }{c}\frac{dn}{d\omega }+\frac{n}{c}.$$
(14)
Using this in the above result for $`\delta `$$`\varphi `$ gives

$$\delta \varphi =\frac{g\omega }{c^2}h(\frac{\omega }{c}\frac{dn}{d\omega }h+\frac{nh}{c}\frac{nh}{c})w\frac{g\omega h^2}{2c^3}\frac{dn}{d\omega }$$
(15)
and

$$\delta \varphi =\frac{g\omega ^2h^2}{2c^3}\frac{dn}{d\omega }.$$
(16)
In principle, this phase shift can be measured relative to light traveling through a fiber with different dispersion to determine or place a limit on $`\alpha `$ in equation (2). Alternatively, the information could be used to make a measurement of g. In fact, a net phase shift can be preserved while the light is returned to the original height by making the return path out of fiber with a different dispersion, as shown in Figure 1. Unfortunately, the phase shift expected per unit length of fiber is small for reasonable values of dispersion, making it necessary to have a relatively long path length to produce a measurable effect. The experimental difficulty is to maintain sufficient light at the end of the path to make a measurement while preserving the phase shift or timing well enough to make the measurement meaningful.

We have considered two experimental models in evaluating the feasibility of a high precision measurement of the gravitational redshift of light using this effect. First we considered using long counter-propagating pulses in a vertical fiber optic Sagnac interferometer. The second model uses optical solitons propagating in two vertically separated fiber optic storage rings. In spite of being less conventional, the latter approach presents fewer technical challenges and is more likely to yield a robust measurement. Since we believe the second of these models is more experimentally feasible, we present it first.

## III Experimental concept: vertically separated optical soliton storage rings

A coherent beam of light pulses is split into two streams, each of which is transformed into a bitstream of solitons. The streams are injected and stored for a period in separate fiber optic storage loops  sitting at different gravitational potentials. The group velocity of the solitons in the two storage rings will differ due to the interplay of dispersion with the gravitationally-induced change in the carrier frequency of the solitons. Unfortunately, the substantial timing jitter that normally accompanies solitons in storage rings makes it impractical to measure directly the relative quantum phase difference between the two beams. Rather, a relative net timing difference between the pulses in the two beams will develop over time, which is proportional to the gravitational frequency shift of the light. This timing difference is used to measure, or place a limit, on $`\alpha `$.

### A Conceptual details

Consider an optical pulse that is split into two coherent pulses with the same intensity profile. Let one travel horizontally for a distance h and the other travel vertically a distance h. The group velocity dispersion parameter $`\beta `$2 for each pulse is given by

$$\beta _2=\frac{1}{2\pi }\frac{dv_g}{d\nu },$$
(17)
where vg is the group velocity of the carrier frequency of the pulse . Assuming the fiber to be identical for each path, there is a gravitationally-induced group velocity difference,

$$\left|\mathrm{\Delta }v_g\right|=\mathrm{\Delta }\nu \left(\frac{dv_g}{d\nu }\right)=2\pi \beta _2\frac{gh\nu }{c^2},$$
(18)
between the pulses as they reach the end of their respective paths, where g is assumed constant over the height difference. (In practice, the variation in g must be measured and taken into account.) At this point, if the pulse streams were each allowed to travel a horizontal distance d along identical fibers and then recombined (in such a manner that the path lengths were identical), there would be a relative pulse timing difference between the pulses . Unfortunately, the effect is quite small for reasonable input parameters and practical distances.

One way of increasing the size of the effect is to increase the distance d. In this model, that is done by using optical soliton pulses and injecting each stream of pulses into ‘identical’ fiber optic storage loops, where amplification is provided by erbium doped fiber amplifiers (EDFA’s). The soliton character prevents dispersion from destroying the pulses, while amplified spontaneous emission noise is controlled using amplitude modulation, and pulse timing jitter is reduced with frequency filtering. This model is schematically illustrated in Figure 2. This and similar techniques for pulse memory and storage and long-haul transmission have been studied extensively in the last decade . A scheme similar to what we are proposing was used in 1993 to store a 10 GHz bit pattern for distances up to 180x106 km (storage times of 15 minutes), as described in reference . The authors of that work conclude, “error free transmission over unlimited distances is possible with this technique.”

The goal of the telecommunications and pulse storage research on optical storage loops is to preserve the information, i.e., a high rate bit pattern, over time and distance. Our goal is somewhat broader, as we want to preserve the bit pattern as well as a small relative overall group velocity difference between the bits in two separate streams. If preserved over time, this relative group velocity difference will be translated into a measurable pulse timing difference. There are sources of random pulse timing jitter, as we will discuss below. However, these noise sources will cause the solitons to random walk in relative timing, while the gravitationally-induced relative velocity difference should provide a small, but inexorable, bias in one direction. There are other effects (as discussed below) which can cause a unidirectional shift in the group velocity of the solitons. In this model, the gravitationally-induced effect is pulled out from the background through differential measurements of the relative pulse timing between configurations with loops at different heights.

Figure 2 shows the basic components of a soliton storage ring experiment. Soliton pulses are used because these solutions to the nonlinear Schroedinger equation balance dispersion against the non-linear Kerr effect, and yield a pulse that is largely immune to dispersive effects in the fiber. Light loss due to scattering is replaced with continuous or lumped amplification by EDFAs. Synchronous amplitude modulation limits the growth of amplified spontaneous noise that would otherwise overwhelm the signal, and a frequency filter is used in concert with the modulator to stabilize the pulse from amplitude and width fluctuations.

### B Noise sources

There are numerous things that can effect the pulse timing in such an optical storage loop.

* Gordon-Haus Jitter
All optical amplifiers introduce noise into the system by amplified spontaneous emission (ASE). In non-linear fiber optics, this ASE noise is also called Gordon-Haus Jitter. When an optical amplifier introduces ASE, the pulse undergoes small amplitude and phase fluctuations at the amplifier output, which lead to a change in the carrier frequency, and therefore the group velocity of the pulse. Over a span of fiber optic cable, this leads to a random jitter in the timing of the soliton pulse. The effect of this jitter on soliton stream transmission has been quantified as

$$BL_T<\left[\frac{\pi f_j^2}{q_0\lambda h\alpha \gamma D}\right]^{1/3},$$
(19)
where B is the bit rate, LT is the length traveled, fj is the fraction of the bit slot that the error should be smaller than, q0 is the initial soliton separation in normalized units, $`\lambda `$ is the wavelength of light, $`\alpha `$ is the attenuation coefficient, $`\gamma `$ is the nonlinearity coefficient, and D is the dispersion coefficient.

Agrawal gives an illustrative example in reference . A system using 1550 nm light with q0 = 10, $`\alpha `$ = 0.2 dB/km, $`\gamma `$ = 10 W-1/km, and D = 2 ps/(nm-km), and assuming fj= 0.2 for tolerable jitter. leads to

$$BL_T<34,000(Gb/s)km.$$
(20)
A 10 Gbit/s system is limited to a transmission length of 3400 km in this example. This length is small compared to the soliton transmission lengths required by our experiment.

Fortunately, there are a number of ways to overcome this bit-length limitation    , including the use of optical filters and external modulators. Filters reduce the overall noise level by eliminating the ASE-noise outside the filter bandwidth and synchronous modulation forces the soliton to realign itself within the bit slot. Both techniques reduce the jitter, and when used together properly, can lead to an unlimited propagation distance.
* Raman Self Frequency Shift (RSFS)
As the soliton pulse travels down a fiber optic cable, the photons scatter off of atoms in the cable. In the process, some of the photon energy can be lost to atomic vibrational modes, yielding an outgoing photon frequency that is lower than that of the incoming photon. This effect, known as the Raman effect, causes a continuous flow of energy from the high frequency part of the soliton pulse to the low frequency part of the pulse. Thus, there is a downshift in the carrier frequency of the soliton as it moves along the cable, which effects the group velocity of the soliton and could give a relative timing difference in our experiment if the effect is asymmetric between the two rings  .

It is well known that a light pulse and its frequency range are related to each other via the Fourier transform. Therefore, the shorter the soliton pulse, the wider its frequency range. Solitons of $``$1 ps or less have a wide enough frequency range that the frequency shift of the soliton due to the Raman effect is very noticeable. Models have shown  that the redshift associated with the soliton scales T0-4, where T0 is the soliton pulse width. So, the noticeable effect drops off quickly for pulse widths above 1 ps. This is an important design consideration in our experiment.
* Soliton Interactions
The same non-linearity that causes solitons to be stable also causes an interaction force between solitons. The interaction force between solitons depends on the separation of the solitons, their relative amplitude, and their relative phase . Only in the case of two solitons with the same relative phase is the soliton-soliton interaction attractive. In all other cases, where the solitons have a different relative phase, the interaction force is repulsive. In the specific case when the relative phase $`\theta `$ = 0 and the relative amplitude r = 1, the relationship between soliton separation and distance traveled down a fiber optic cable is given by

$$\mathrm{exp}[2(qq_0)]=\frac{1}{2}[1+\mathrm{cos}(4\xi e^{q_0})],$$
(21)
where q0 is the initial soliton separation, and q is the soliton separation at a distance $`\xi `$ along the fiber. From this equation it can be seen that the solitons periodically combine and then separate once again with a period $`\xi `$p, given by

$$\xi _p=\frac{\pi \mathrm{sinh}(2q_0)\mathrm{cosh}(2q_0)}{2q_0+\mathrm{sinh}(2q_0)}.$$
(22)
It is important to note the interaction force drops exponentially with the distance between solitons.

Soliton interactions for two solitons, three solitons, and a random sequence of solitons have been studied numerically . It was found that in a given soliton stream with a number of consecutive bits filled, it only the first and last solitons of the consecutive bit stream are effected by the soliton interaction because the interior solitons experience approximately equal forces from each side.
* Third Order Dispersion
Dispersion in fibers is usually described by taking the lowest order term in a group velocity dispersion (GVD) expansion. In most cases, this term is dominant and the first order expansion is a good approximation. However, near the zero point of dispersion, the next term in the expansion, known as third order dispersion can be of comparable magnitude to the lowest order term, and must be considered . This effect has attracted attention recently because it can be a significant effect in a dispersion-managed system, where the fiber alternates between being dispersive and anomalously dispersive and the average dispersion is slightly anomalously dispersive and very small . In general, this effect is small unless the soliton frequency is very near the zero of dispersion ($`\beta `$2) for the fiber.
* Polarization mode cross coupling
Even in a single mode optical fiber, two polarization modes can be maintained. Dispersion between these two modes is known as Polarization Mode Dispersion (PMD.) The internal structure of a soliton is very stable under the effects of PMD. However, the different polarization modes most often have different linear and non-linear indexes of refraction and therefore travel at different group velocities. Over a long span of fiber, even a pulse that starts initially in one polarization mode may be split into both modes via PMD and therefore reach the end of the fiber at different times. In addition, noise involved in PMD causes a soliton timing jitter  .
* Temperature fluctuations Temperature fluctuations can cause a substantial change in the optical path of a pulse during storage in a loop. This effect can be large compared to the delay induced by the gravitational frequency shift and large compared to the other effects listed above.

### C Design issues and performance

Careful design can help control the various sources of noise listed above. The Raman self frequency shift and soliton-soliton interactions can be limited by using relatively wide solitons (\>1 ps) at a relatively low bit rate (10 Gbit/s). Also, soliton-soliton interactions can be reduced by the use of synchronous amplitude modulation  and perhaps amplitude variation. In addition, one can reduce the effect of soliton-soliton interactions by using the relative timing between interior bits only . Third order dispersion can be made negligible by staying away from the zero of dispersion. In fact, the signal in this model depends on the solitons seeing a net dispersion as they transit the loop. Polarization mode cross coupling can be reduced through the use of polarization maintaining (PM) fiber in the loop. However, it would increase the fiber cost substantially if a large fraction of the loop were constructed with PM fiber. Some of the effect of polarization mode cross coupling can be removed by making a differential measurement. In this model it would be necessary to place polarization controllers in strategic locations throughout the circuit. In addition, it might be possible to use polarization specific filtering and detection or make use of trapped or vector solitons (where the two polarization modes have the same group velocity). The timing clock signal for the synchronous modulation can be extracted from the transmitted solitons themselves, following Nazakawa et. al.

Even taking into consideration proper design, there will exist residual relative timing jitter and shifts from the noise sources above. We can place a reasonable limit on the size of the random jitter by noting the results of previous bit storage experiments. Nazakawa et al. established long-term storage of 24 ps wide solitons at 10 Gbit/s . Numerous other groups have transmitted soliton bit patterns for shorter distances at higher (20-160 Gbit/s) rates . Liao and Agrawal have shown through numerical simulations that it should be possible to limit the jitter to approximately 1 ps at a rate of 40 Gbit/s . Based on these results, it seems reasonable to assume the random jitter can be controlled to better than the 10 ps level. This corresponds to an uncertainty in the pulse position at detection of approximately 2 mm, ignoring thermal effects for the moment.

Now, let us return to our experimental model and consider only noise that can cause random jitter. A 10 Gbit/s soliton beam propagating in a fiber with D=+15 ps/nm-km ($`\beta _2=19ps^2/km`$) is split into two streams. One travels vertically a distance h and is injected into a fiber optic storage ring, while the other is injected into a similar fiber optic storage ring at the height of the original soliton pulse source. There will be a gravitationally-induced shift in the group velocity of the solitons that move vertically, according to equation (9). If the storage rings are constructed of similar fiber (assume D=+15 ps/nm-km, SMF-28 fiber), the average, relative group velocity difference between the streams of solitons traveling in the two storage rings will remain. Over a storage time T, the light will travel a distance d in the storage loops and the group velocity difference will induce a relative shift in the length traveled along the two paths (in meters) of

$$\delta d=1\times 10^{10}hT.$$
(23)
There will be

$$N=\frac{Ln\left(rate\right)}{c}=50L$$
(24)
solitons in each storage ring of length L, assuming the rings are filled. The statistical average error on the relative position of solitons in one ring as compared to the other when the beams are recombined is

$$\delta s=\frac{0.002}{\sqrt{N}}meters.$$
(25)
Figure 3 shows $`\delta `$d as a function of vertical loop separation d, and storage time T. The sensitivity to $`\alpha `$, given by $`\delta `$s/$`\delta `$d, is plotted in Figure 4. This figure shows an interesting limit on $`\alpha `$ can be achieved under circumstances where only random jitter noise effects are considered.

Unfortunately, unidirectional soliton timing shifts within each ring are unavoidable and quite large. However, so long as these shifts are reproducible and consistent for a given ring, it should be possible to extract the gravitationally-induced timing shift using differential measurements between different loop heights. Our conceptual design creates a stable environment for the storage loops in hopes of achieving this reproducibility.

Figure 5 illustrates the basic concept of making differential measurements in this experimental model. The fiber storage rings are constructed as “modules” that can be moved carefully from place to place without disturbing the active elements of the storage loop or the environment that surrounds them. Fiber leads to and from each module are detachable, so the relative spatial configuration of the two modules can be changed in order to make differential measurements. Between separate differential measurements there will be relative changes in the input and output fibers. However, the light spends the vast majority of the storage time inside the loop. The input and output fibers constitute an extremely small fraction of the total path length, and will have a correspondingly small effect on the pulse timing.

Figure 6 shows the basic layout of the experimental optical circuit in more detail. The source and both storage rings (and detectors, if required) are in stable environments. The ring modules are constructed to identical, and can be disconnected and moved to different heights, as needed.

Systematic effects, such as those arising from small differences between the modules/storage rings, can be taken out with differential measurements. Each pass around the ring, a small fraction of the power is picked off so that the relative timing in the two rings can be compared.

Figure 7 shows one concept of a fiber optic storage ring module. It consists of the basic storage ring components (fiber, bandpass filter, EDFAs, and amplitude modulator) in a thermally isolated compartment that can be moved as needed. The fiber is coiled on spools that have low thermal expansion coefficients. The coils are arranged so the net area subtended by the loop is very small. This will reduce any effects due to the rotation of the earth during the measurement.

The module is placed in a carefully controlled, stable thermal environment. The temperature inside the module is carefully monitored. If necessary, this information can be used for a feedback system to keep the temperature stable over time. The absolute accuracy of the temperature sensors is not critical. However, they need to give reliable relative measurements over time.

Similarly, the pressure in each module is monitored and adjusted slightly as needed. The relative pressure between the inside of the module and the outside is relatively small in order to minimize difficulties in module enclosure construction. We plan to seal the modules while at the same height and pressure and then transport them to different vertical positions while maintaining a constant internal pressure.

Unfortunately, relatively small fluctuations in the ambient environment of the loops can cause large effects over such long storage times. For example, assuming a coefficient of thermal expansion for the fiber of approximately 10-7/C and a relative temperature difference between two storage loops of $`\delta `$t C leads to a fractional optical path length difference of 10$`{}_{}{}^{7}\delta `$t. This compares to the gravitationally induced effective fractional path length difference of approximately 10-19h, for a loop height separation of h meters. It is impractical to assume the temperature in the modules can be controlled to high enough accuracy to make these numbers even remotely comparable. Rather, it will be necessary to calibrate out effects due to temperature and pressure variations. One way of doing this is to use wavelength division multiplexing in the soliton bit stream. For example, suppose the soliton bit stream contains pulses centered at the zero-dispersion frequency of the fiber as well as pulses in the dispersive region. The relative pulse timing between the zero-dispersion-centered solitons gives a measure of the relative environmental fluctuations since the gravitational effect is nonexistent at the zero dispersion frequency. The relative pulse timing between the pulses in the dispersive region will contain the environmental fluctuations and the gravitational effect. In principle, the gravitational effect can be isolated by subtracting out the relative timing difference of the zero-dispersion-centered solitons. The ultimate success of this experiment hinges, in part, on how well this can be done.

All wires and fibers leading into the module are detachable, so the modules can be exchanged with ease, without disturbing the active elements of the loop. The items inside the module are secured so that minor movements of the module leave no net change in the positions or stresses inside the module.

A careful survey is done to locate the modules with respect to each other and to positions where measurements of the local gravitational field strength are done using a reasonably accurate gravimeter (good to 1 part in 106 or 107).

Each module rests on a stable optical table. The relative position of the two modules is monitored during a store via an interferometric system (or repeated optical surveys if such a system is simpler and more cost effective). Significant changes, perhaps due to diurnal variation in the building, can be taken out via small, remote control, height adjusters located underneath one of the modules. If necessary, these adjustments could be made automatically with feedback from the interferometric system. Accuracy on the order of 200 microns should be achievable and sufficient.

The change in the group velocity, and therefore the signal size, is directly proportional to the group velocity dispersion in the fiber loop. It is necessary to measure this dispersion as precisely as possible so it does not limit the measurement precision on $`\alpha `$. Measurements of fiber chromatic dispersion to $`\pm `$0.02 ps/nm-km have been reported . Fujise et al. use a wavelength-division-multiplexing (WDM) phase shift method which compares the relative phase shift between modulated light of different wavelengths passing through the fiber simultaneously. They achieve a measurement stability of $`\pm `$0.02 ps/nm-km over a ten-hour period. We will need to improve on this precision, or develop the ability to correct for changes over time using this (or a similar) technique in situ as a monitor. Another method we are considering is to send separate soliton pulses with different carrier frequencies through the loop and look at the relative timing between the pulses before and after the loop. The stability of our dispersion measurement is intimately tied in with the stability of the overall timing at the end of a store and our ability to extract the gravitational effect using differential measurements.

### D Potential for a flight project

Since the signal size of the soliton ring model depends linearly on the height difference of the two loops in a given gravitational field, it is desirable to make this quantity as large as possible. To this end, it is interesting to explore the potential gains and difficulties that might be encountered performing the soliton storage ring experiment in orbit on the space shuttle, or the International Space Station. In such an experiment, one storage loop would remain on the shuttle/station, while the other loop would be transported to an orbit with a different radius or height, but the same angular velocity. In this scheme, loop height differences of tens of kilometers or more could be achieved.

The acceleration due to gravity in orbit will differ from that on the surface of the earth by about 10%. (The value of g is around 8.94 m/s2 at an altitude of 160 nautical miles, which is the altitude of STS-101.) This comes in the signal with a linear dependence and must be taken into account for an exact calculation. However, for simplicity, we will ignore this difference, along with the change in g as a function of height along the loop. In practice, these (O(10%) effects) must be measured and taken into account. With this simplification, equation (14) is still valid. Figure 8 shows the potential sensitivity to $`\alpha `$ for heights of 10, 20, and 50 km, respectively.

Conceptually, one module remains on the shuttle (or space station), while the other is ensconced in a small, unmanned vehicle that contains environmental control, power for EDFA pump lasers and thermal control, environmental monitoring, devices capable of relating the exact position of the one module with respect to the other, and minor propulsion and navigational capability for stability (and possible maneuvering, though perhaps this is better done by the shuttle). The two modules must be tethered together by the fibers that carry the soliton stream to and from the second module. The fibers should be surrounded by a thin, flexible covering that evens out the azimuthal asymmetry in the thermal environment and provides some degree of stress relief.

Conceivably, the solitons could be generated separately on the two modules with the timing communicated via light pulses sent back and forth. This has the advantage of removing the need for the fiber tether and allowing the possibility of a greater height difference, and perhaps a sensitivity to higher-order terms. Unfortunately, this raises a host of other issues, such as relative differences in the loop soliton sources and gravitational effects on the clock rate of electronics that communicates between the modules.

It is also important to note that the primary elements of this experiment are small, light, and mechanically robust. In addition, the experiment requires a relatively modest amount of power. As such, it could be flown on a variety of platforms.

This experimental model is quite novel. It would have very different systematic errors than the clock-in-space experiments that are underway.

## IV Experimental model of vertical fiber optic Sagnac interferometer with varying dispersion

Coherent light traveling from a source is split into two beams that travel in opposite directions around a loop or coil of optical fiber. After traversing the coil, the two beams of light are recombined and the phase difference between the clockwise (CW) and counter-clockwise (CCW) beams is measured. This technique is used widely to make high precision fiber optic gyroscopes (FOG) that make use of the fact that a phase difference between the two paths is induced by a rotation of the coil . In our model, the coil lies in a vertical plane. Unlike the FOG, the fiber dispersion must vary asymmetrically around the circumference of the loop in order yield a net phase shift due to the interplay of dispersion and the gravitational frequency shift. At the top and bottom of each loop in the coil are splices to join the fibers of differing dispersion. The pattern of varying dispersion is such that the CW (CCW) beam always travels upward through fiber with positive dispersion and downward through fiber with negative or anomalous dispersion. The situation is reversed for the CCW (CW) beam.

As the CW and CCW beams propagate through the coil, a phase difference between the beams develops with time due to the gravitational frequency shift. This phase difference is given by

$$\delta \phi =\frac{2\pi Ngh^2\nu ^2}{nc^3}\left[\left(\frac{dn}{d\nu }\right)_L\left(\frac{dn}{d\nu }\right)_R\right]$$
(26)
for beams of initial frequency $`\nu `$, traveling N turns around a loop of height h in a uniform gravitational field with strength g. The dispersion on one side of the fiber loop is (dn/d$`\nu `$)L while that on the other side is (dn/d$`\nu `$)R.

At 1550 nm, for an interferometer constructed of fibers with dispersion +15 ps/nm-km (=3.7x10-17 s) and
-90 ps/nm-km (=-2.2x10-16 s), respectively, and setting n=1.5, equation (26) becomes

$$\delta \phi =\left(1.6x10^{11}\right)Nh^2.$$
(27)
Figure 9 shows the size of the phase shift as a function of total path length and loop height. Although many sources of noise need to be considered, the measurement precision is limited ultimately by light loss in the fiber and at the splices.

In order to estimate, roughly, the potential sensitivity of such an experiment, it is necessary to fold in the shot noise. The photon counting noise in the phase can be estimated by the Heisenberg uncertainty principle as

$$\left(\delta \phi \right)_{shot}\frac{1}{\sqrt{N_{photons}}}.$$
(28)
For a 10 mW laser source operating at 1550 nm, this estimate becomes

$$\left(\delta \phi \right)_{shot}=\frac{1}{\sqrt{\left(5x10^{17}\right)te^{\alpha d}}},$$
(29)
where t is the signal integration time (in s), $`\alpha `$ is the total fiber light loss (in dB/km), and d is the total distance traveled by the light (in km). Assuming a fiber light loss of 0.5 dB/km and a loss at each splice of 0.5 dB yields the shot noise limit shown in Figure 9.

The results in figure 9 are from analytical calculations. We have also used the split step Fourier method to solve the non-linear Schrodinger equation for this model numerically . The results, shown in Figure 10, are very similar to those obtained by analytical calculations. The difference in the shot noise curve comes about because the analytical model treats the light loss at splices as a continuous process, while the numerical model incorporates discrete losses at the splices and should be more accurate.

Things besides shot noise can limit the sensitivity of this experiment. Source relative intensity noise, thermal fluctuations, changing mechanical stresses, polarization mode cross coupling, and coherent backscattering at the splices can lead to nonreciprocal phase shifts. In addition, electronics noise can be an issue. Without constructing a very detailed model and/or building prototypes, it is difficult to estimate the level of these sources of noise in this experimental model. One useful point of comparison is the fiber optic gyro, mentioned earlier. The FOG is a continuous wave fiber Sagnac interferometer that is vulnerable to many of the same sources of noise. Very high sensitivity FOGs are sensitive to a rotation-induced phase change of 1x10-9 radians. A sensitivity of one $`\mu `$rad is considered difficult, but achievable  .

Figure 11 shows the signal-to-noise ratio for this model with varying loop height, shot noise and a 1 $`\mu `$rad noise floor. Without using squeezed light, the shot noise provides a hard cutoff in the path length for a given loop height. An optimistic estimate of what could be done using squeezed light is shown in Figure 12. For this curve, the shot noise is reduced by a factor of 100 .

In our view, Figures 11 and 12 represent an optimistic assessment of what might be achievable with current technology. Our model differs from a FOG in some very significant ways. FOG designers struggle to maintain reciprocity between the counter propagating beams in order to achieve the high sensitivity of the best instruments. Non-reciprocal dispersion is a necessary component of this model. In addition, FOGs are wound about titanium cores in a carefully maintained and monitored environment. In our model, the loop necessarily extends a large vertical distance, which would complicate attempts at noise reduction.

Another weakness of this model is that it would be very difficult to make a differential measurement, which would be helpful for extracting the signal. It is hard to imagine being able to rotate the device or change the height without causing arbitrary phase changes between the counter-propagating beams. The growth of the phase shift with time, compared in different orientations, would be the best way to pull out the signal. This is another place where our model differs significantly from the FOG. FOGs work with continuous wave light and our model requires long pulses. 0

To summarize our work on this model, it is conceivable that a high precision measurement of the gravitational redshift of light could be made using a vertical fiber Sagnac interferometer with varying dispersion. However, the technical difficulties are challenging and further modeling and prototyping would be useful in evaluating the eventual sensitivity that might be achieved.

## V CONCLUSIONS

The interplay of optical dispersion and the frequency shift of light due to a changing gravitational potential can be used to test the equivalence principle for gravitation by converting the frequency shift into a phase shift or time delay. It is possible an experiment using solitons in optical fiber storage loops could be performed that uses this effect to make a high precision measurement of the gravitational redshift of light.

## VI Acknowledgements

We are indebted to the following people for useful discussions: Govind Agrawal, Robert Boyd, Jeffrey Clark, Turan Erdogan, Alan Evans, John Heebner, Daniel Kerr, Jeff Korn, Taras Lakoba, Adrian Melissinos, John Moores, Michael Perlmutter, Sarada Rajeev, and Ram Yahalom. This work was supported by the University of Rochester.

## VII References


format: text

text:
# Consequences of a 𝑍₂ Symmetry for Neutrino Oscillations

## Abstract

A $`Z_{2L}\times Z_{2R}`$ generation symmetry in the neutrino sector predicts the atmospheric neutrino mixing to be maximal, and the MNS matrix element $`U_{e3}`$ to be zero, consistent with observations. Solar neutrino mixing may be maximal but is not required by the symmetry. Neutrino masses of the first two generations are predicted to vanish, providing a first approximation to the oscillation data. The consequence of a smaller $`Z_2`$ symmetry is also discussed. In that case, deviation from the $`Z_{2L}\times Z_{2R}`$ result is of the order of the neutrino mass ratio between the first two generations and the third generation.

Right-handed neutrinos must be present to explain the deficit of solar and atmospheric neutrinos through oscillations. These are unique probes at high energy for they have no standard-model quantum numbers to allow a pollution by standard-model interactions. Recent results from Super-Kamiokande suggest that the mixing for solar neutrinos and the atmospheric neutrinos are both maximal, and the $`(\mathrm{\Delta }m)^2`$ for solar neutrino oscillation is much smaller than the $`(\mathrm{\Delta }m)^2`$ for atmospheric neutrino oscillations . The negative result of CHOOZ  also limits the magnitude of the MNS mixing matrix  element $`|U_{e3}|^2`$ to be smaller than 0.02 to 0.047, depending on the exact value of the atomospheric neutrino mass difference. Is there a simple way to understand these facts? Many models have been proposed . We suggest that a $`Z_{2L}\times Z_{2R}`$ generation symmetry in the sea-saw scenario will naturally explain many of these observed results. By $`Z_2`$ we mean a finite group with elements $`\{e,g\}`$ so that $`g^2`$ is the identity element $`e`$. Specifically, under such a symmetry, the atmospheric neutrino mixing is maximal, and the MNS matrix element $`U_{13}U_{e3}`$ is zero. The neutrino masses of the first two generations vanish, which may be a fair approximation to Nature since the mass difference observed in atmospheric neutrino oscillation is much larger than that seen in solar neutrino oscillation. Under such a symmetry, the solar neutrino mixing may be maximal but is not required to be so.

Let $`D`$ be the Dirac mass matrix for neutrinos in the flavor basis, and $`M`$ the Majorana mass matrix in the same basis for right-handed neutrinos. We shall assume oscillations to occur only among the active neutrinos, so both are $`3\times 3`$ complex matrices, although $`M`$ is symmetrical. When the right-handed neutrinos are integrated out, the effective mass matrix for the left-handed neutrinos in the flavor basis is

$`m^{}=D^TM^1D.`$ (1)
To diagonalize it, we must use the MNS matrix $`U`$ to rotate the flavor basis into the energy eigenbasis, then

$`m^{}=U^{}mU^{},`$ (2)
with $`m=\mathrm{diag}(m_1,m_2,m_3)`$.

Let $`\nu _{aL}`$ and $`\nu _{aR}`$ be respectively the left- and right-handed neutrino field of the $`a`$th generation. Assume the mass terms in the neutrino Lagrangian to be invariant under the $`Z_{2L}\times Z_{2R}`$ transformation $`\nu _{2A}\nu _{3A},\nu _{3A}\nu _{2A}`$, and $`\nu _{1A}\nu _{1A}`$, where $`A`$ is either $`L`$ or $`R`$. Then the matrices $`D`$ and $`M^1`$ must be of the form

$`D`$ $`=`$ $`m_D\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 1\\ 0& 1& 1\end{array}\right),`$ (3)
$`M^1`$ $`=`$ $`m_M^1\left(\begin{array}{ccc}\mu _1& \mu _2& \mu _2\\ \mu _2& 1& \mu _3\\ \mu _2& \mu _3& 1\end{array}\right).`$ (5)
Using (1), this leads to

$`m^{}={\displaystyle \frac{2m_{D}^{}{}_{}{}^{2}(1\mu _3)}{m_M}}\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 1\\ 0& 1& 1\end{array}\right).`$ (6)
Other than the common scale factor which may be complex, the matrix is real and symmetrical, so it can be diagonalized by a real orthogonal transformation. The eigenvalues are $`(0,0,1)`$, but the matrix $`U`$ used for diagonalization is not unique because of degeneracy of the first two eigenvalues. For the purpose of later generalization it is useful to deduce these results explicitly from

$`(m)_{ij}`$ $`=`$ $`U_{ai}(m^{})_{ab}U_{bj}`$ (7)
$`=`$ $`{\displaystyle \frac{2m_{D}^{}{}_{}{}^{2}(1\mu _3)}{m_M}}\left(U_{2i}U_{3i}\right)\left(U_{2j}U_{3j}\right).`$ (8)
In order for $`m`$ to be diagonal, two of the following three combinations must vanish: $`U_{21}U_{31},U_{22}U_{32}`$, and $`U_{23}U_{33}`$. Let us first consider the case when $`U_{23}U_{33}0`$, but the other two combinations vanish. In that case $`m_1=m_2=0`$ and $`m_3=(2m_{D}^{}{}_{}{}^{2}(1\mu _3)/m_M)(U_{23}U_{33})^2`$, thus giving rise to the extreme limit of a normal hierarchy for the neutrino masses, which is a good approximation to reality because $`(\mathrm{\Delta }m)_{12}^2(\mathrm{\Delta }m)_{23}^2`$. To simplify writing let $`a=U_{21}=U_{31}`$, $`b=U_{22}=U_{32}`$. Then $`U`$ is of the form

$`U=\left(\begin{array}{ccc}c& d& e\\ a& b& f\\ a& b& g\end{array}\right).`$ (9)
We will adopt the usual phase convention so that $`c,d,f`$ and $`g`$ are real, and the imaginary part of $`a`$ and of $`b`$ are proportional to that of $`e`$. Normalization of the second and third rows of the matrix $`U`$ requires $`|a|^2+|b|^2+f^2=1=|a|^2+|b|^2+g^2`$. Hence $`f`$ and $`g`$ have the same magnitude. Orthogonality of the second and the third rows of $`U`$ then requires $`f=g`$ and $`f^2=|a|^2+|b|^2=\frac{1}{2}`$. In other words, atmospheric neutrino mixing is maximal, consistent with the Super-Kamiokande observation. Using the fact that the first row of $`U`$ must be orthogonal to both the second and the third roles, we conclude that $`e=0`$, consistent with the CHOOZ observation. In that case there will be no CP violation observable through neutrino oscillations. With the present phase convention both $`a`$ and $`b`$ become real, and the magnitudes of $`a,b,c,d`$ are related by unitarity. So, there is only one free parameter $`\theta _{12}`$ left in describing $`U`$. For $`c=d=1/\sqrt{2}`$, maximal mixing occurs in solar neutrino mixing, but this is not required to be so by the symmetry. In summary, judging from neutrino oscillations, the $`Z_{2L}\times Z_{2R}`$ symmetry appears to be a good approximation to Nature.

The other two solutions of (8) leads to either $`m_10`$ or $`m_20`$, with the other two diagonal matrix elements of $`m`$ zero. The former case leads to $`U_{11}=0`$, and the latter case leads to $`U_{12}=0`$, neither is consistent with solar neutrino and CHOOZ observations. They will therefore be rejected.

Note also that the symmetry interchanges generations 2 and 3, with a minus sign. If instead it interchanged generations 1 and 2, which at first sight might seem a more reasonable thing to do in view of the mass hierarchy, the result would be the same as above but with the first and third rows of the matrix $`U`$ permuted. This would not be consistent with experiment.

If the charged lepton and the quark Dirac mass matrices are subject to the same symmetry, then all of them would have a mass spectrum proportional to (0,0,1), not a bad first approximation. In that case the CKM matrix would be $`\mathrm{𝟏}`$, again a reasonable first approximation given that the Cabibbo angle is small.

There is the question of how a small but non-zero value of $`m_{1,2}`$ can be obtained. One possibility is to have a small breaking of $`Z_{2L}\times Z_{2R}`$, into a diagonal $`Z_2`$ where the left-handed and the right-handed particles are simultaneously transformed. This does not alter the parametrization of $`M`$, but allows small parameters to occur in $`D`$ where its present matrix elements are zero. It should be possible to tune these parameters to get finite masses for the first two generations. Since a different set of parameters may be present for the charged leptons, the up-type quarks, and the down-type quarks, it is quite conceivable that different mass spectra can be obtained for these different particles, while obtaining also a realistic CKM matrix.

Let us now explore this scenario more systematically. The $`Z_2`$ symmetry requires the (symmetric) matrix $`m^{}`$ in (1) and (2) to obey $`m_{12}^{}=m_{13}^{}`$ and $`m_{22}^{}=m_{33}^{}`$. Using (2), these two requirements can be written as

$`{\displaystyle \underset{i}{}}\alpha _i\beta _im_i`$ $`=`$ $`0,`$ (10)
$`{\displaystyle \underset{i}{}}\gamma _i\beta _im_i`$ $`=`$ $`0,`$ (11)
where

$`\alpha _i`$ $`=`$ $`U_{1i}^{},`$ (12)
$`\beta _i`$ $`=`$ $`U_{2i}^{}U_{3i}^{},`$ (13)
$`\gamma _i`$ $`=`$ $`U_{2i}^{}+U_{3i}^{},`$ (14)
thus giving rise to constraints on $`U`$ controlled by the neutrino mass values. For the extreme hierarchy spectrum $`(m_1,m_2,m_3)(0,0,1)s_1`$ discussed above, the constraints are $`\alpha _3\beta _3=\gamma _3\beta _3=0`$. Being second order equations there are two solutions: (1a) $`\alpha _3=\gamma _3=0`$, and (1b) $`\beta _3=0`$. Solution (1a) is the same as the $`Z_{2L}\times Z_{2R}`$ solution before. To see that, note that the third column of $`U`$ is fixed by this constraint and unitarity to be $`(U_{13},U_{23},U_{33})=(0,f,f)`$, with $`f^2=\frac{1}{2}`$. Orthogonality of the first two columns to the third then fixes $`U`$ to be of the form (9), with parameters identical to those in $`Z_{2L}\times Z_{2R}`$ as per the arguments given below (9). Solution (1b) gives $`U_{23}=U_{33}`$, which is not very restrictive. For a realistic hierarchical spectrum $`s_2=(\eta ,ϵ,1)`$, with $`ϵ1`$ and $`\eta `$ either of the same order or much less than $`ϵ`$, the constraints become $`\alpha _3\beta _3=O(ϵ)`$ and $`\gamma _3\beta _3=O(ϵ)`$. The two corresponding solutions are: (2a) $`\alpha _3=O(ϵ)`$, $`\gamma _3=O(ϵ)`$, and $`\beta _3=O(1)`$, and (2b) $`\beta _3=O(ϵ)`$ and $`\alpha _3,\gamma _3=O(1)`$. As before, solution (2b) is not very restrictive. For (2a), $`|U_{e3}|`$ is of order $`ϵ`$ and the atmospheric neutrino oscillation is close to maximal, with deviation of the order of $`ϵ=m_2/m_3`$. In this way a $`Z_2`$ symmetry with a realistic mass spectrum as an input gives rise to practically the same predictions as the larger $`Z_{2L}\times Z_{2R}`$ symmetry.

Whether a more realistic mass spectrum should be obtained by an explicit breaking of $`Z_{2L}\times Z_{2R}`$ like the one described above, by renormalization-group corrections, or by something else, will be left to future investigations.

This research is supported in part by the Natural Sciences and Engineering Research Council of Canada, and the Fonds pour la formation de Chercheurs et l’Aide à la Recherche of Québec. I am indebted to T.K. Kuo, Greg Mahlon, Gary Shiu, and Tony Zee for stimulating discussions.


text:
# 85-GHz BIMA observations of the double-hotspot radio galaxy 3C 20

## 1 Introduction

In the beam model for powerful extragalactic double radio sources, hotspots (the bright compact regions at the ends of the source) are the visible manifestation of a strong shock as the relativistic beam of energetic particles is suddenly decelerated by interaction with the external medium surrounding the radio lobes.

This model is challenged by the observation that the lobes of radio galaxies and quasars very frequently have more than one hotspot. To account for this, we must either assume that the beam end-point moves about from place to place in the lobe (the ‘dentist’s drill’ model of Scheuer 1982) or that material flows out from the initial impact point of the beam to impact elsewhere on the lobe edge \[the ‘splatter-spot’ model of Williams & Gull (1985) or the jet-deflection model of Lonsdale & Barthel 1986\]. Both these models predict that one of the hotspots, the one associated with the first termination of the jet, should be more compact than the other or others; it is in fact observed that where jets are explicitly seen to terminate, they always do so in the most compact, ‘primary’ hotspot (Laing 1989, Leahy et al. 1997, Hardcastle et al. 1997). But from single-wavelength radio observations it is very difficult to distinguish between different models of multiple hot spot formation.

One area in which different models do make different predictions is that of the high-frequency spectra of the secondary (less compact) hotspots. If these are relics left behind by the motion of the jet, as in the ‘dentist’s drill’ model, then in general we expect shock-driven particle acceleration to have ceased \[though, as pointed out by Cox, Gull & Scheuer (1991), it is possible for a disconnected hotspot to continue to be fed for some time if the disconnection occurs at a significant distance upstream\]. Synchrotron losses will then deplete the high-energy electrons in the secondary hotspot, so that the secondary hotspot is expected to have a steeper high-frequency spectrum than the primary. If secondary hotspots continue to be fed by outflow from the primary hotspot, then there is still an energy supply and particle acceleration will continue to operate.

To have the best chance of distinguishing between models using these spectral differences we must observe at the highest available frequencies, since the synchrotron lifetime of an electron emitting at a frequency $`\nu `$ is proportional to $`\nu ^{1/2}`$. Observations show that (primary) hotspots often have a spectral cutoff around the mm-wave region of the spectrum, perhaps corresponding to a maximum electron energy attained by the acceleration process. Mm-wave observations therefore provide one of the best opportunities to test the different models of multiple hotspot formation.

Only a few multiple-hotspot sources have been observed at millimetre wavelengths with resolution sufficient to distinguish between the components of the hotspot complex. Two such sources are the nearby luminous classical double Cygnus A, 3C 405 (Wright & Birkinshaw 1984, Wright & Sault 1993) and the peculiar radio galaxy 3C 123 (Looney & Hardcastle 2000), both of which show close hotspot pairs (separated by a projected distance of about 8 kpc and 5 kpc respectively). In both these cases the radio-to-mm spectra of the two hotspots were found to be identical within the errors; there was no evidence for a steeper spectrum in the secondary hotspot.

In this paper we present 85-GHz observations of a well-known double-hotspot radio galaxy, 3C 20. The double hotspot in the E lobe of this $`z=0.174`$ object was one of the earliest to be resolved in radio observations (Jenkins, Pooley & Riley 1977; Laing 1981), because of the components’ large angular separation (6.5 arcsec, corresponding to 27 kpc). It has since been well observed with the NRAO Very Large Array (VLA) at various radio frequencies (e.g. Hiltner et al. 1994, hereafter H94; Hardcastle et al. 1997; Fernini et al. 1997). Its 8.5-GHz radio hotspot structure is shown in Fig. 1. The western hotspot, which lies on what is believed to be the jet side (Hardcastle et al. 1997) and is actually a more compact hotspot pair (or possibly a single hotspot with an extended tail) is known to be a mm-wave and optical source (H94, Meisenheimer et al. 1989, Meisenheimer, Yates & Röser 1997) with a spectrum which cuts off in the optical. But the eastern hotspots were not detected in the optical (H94, Meisenheimer et al. 1997) and until now they have not been observed at high radio frequencies.

Throughout the paper we use a cosmology with $`H_0=50`$ km s-1 Mpc-1 and $`q_0=0`$. With this cosmology, 1 arcsecond at the distance of 3C 20 corresponds to 3.99 kpc. Spectral index $`\alpha `$ is defined in the sense $`S\nu ^\alpha `$, and $`\alpha _{\nu _1}^{\nu _2}`$ denotes the two-point spectral index between frequencies of $`\nu _1`$ and $`\nu _2`$.

## 2 Observations and data analysis

3C 20 was observed in the B and C configurations of the 10-element BIMA array111The BIMA array is operated by the Berkeley Illinois Maryland Association under funding from the U.S. National Science Foundation.. The data were taken on 1999 Oct 11 (C array) and 1999 Oct 30 (B array), with the correlator configured to give two 800-MHz bands centred on 83.15 and 86.60 GHz. These were combined in the final images to give an effective frequency of 84.87 GHz. Combining the B- and C-configuration data gives good coverage of the uv plane between 2.2 and 66 k$`\lambda `$, so the resulting images have a resolution of about 3 arcsec and are sensitive to structures up to $`50`$ arcsec in size, well matched to the largest angular size of 3C 20.

The data were processed using the MIRIAD package (Sault, Teuben, & Wright 1995). Amplitudes were calibrated using observations of Uranus (B-array) and Mars (C-array) to bootstrap the flux density of the phase calibrator, the quasar 0136+478, which was set to 3.68 Jy. The uncertainty in the amplitude calibration is estimated to be about 7 per cent.

Final images were produced using the AIPS task IMAGR, and are shown in Fig. 2. The western hotspot and the two eastern hotspots are clearly detected, and (as shown in Fig. 2b) there is also some low-surface-brightness emission from both the lobes of 3C 20. The weak radio core is not detected at 85 GHz. The net flux density of the source at 85 GHz is about 0.21 Jy. This is somewhat lower than the 90-GHz flux density of 0.51 Jy quoted by Steppe et al. (1988), but the error on the value of Steppe et al., though not given explicitly in their paper, is certainly large.

## 3 VLA observations

To compare our data with observations at longer wavelengths, we obtained existing VLA data or images at 1.4, 4.9, 8.5 and 15 GHz. The 1.4 GHz image was taken from Leahy, Bridle & Strom (1998) based on unpublished observation by R.A. Laing with the VLA A configuration. The 8.5-GHz data were described by Hardcastle et al. (1997) and used the VLA A, B and C configurations. We retrieved 4.9- and 15-GHz snapshot observations taken in 1987-88 from the VLA archive, using B and C configurations and C and D configurations respectively, and calibrated and reduced them in the standard manner within AIPS, using 3C 286 as the primary flux calibrator. The flux calibration is estimated to be accurate to within 2 per cent. \[In the case of the D-configuration 15-GHz observations no scan on a suitable calibration source was made, and we calibrated the amplitude by reference to the radio galaxy 3C 123, whose 15-GHz total flux density we know (Looney & Hardcastle 2000).\] All these datasets have shortest baselines $`2`$ k$`\lambda `$, well matched to the BIMA observations. For comparison with Fig. 2, in Fig. 3 we show VLA images at these four frequencies with a resolution of 3 arcsec.

## 4 Spectra

The low resolution of the BIMA data compared to the size of the compact hotspots of 3C 20 (Fig. 1) means that there is not a single obviously correct procedure for using these data to constrain the hotspot spectra. The 3-arcsecond beam certainly contains more than one spectral component in the case of the compact hotspot in the eastern lobe (hereafter the ‘NE hotspot’) and the compact pair of hotspots in the western lobe (the ‘W hotspot’), because these two components are surrounded by more or less diffuse emission with a steeper spectrum at cm wavelengths (Fig. 4) which may or may not contribute to the 85-GHz flux density. On the other hand, the secondary hotspot in the E lobe (the ‘SE hotspot’) might well be spectrally homogeneous, given its relaxed appearance at high resolution.

We have therefore chosen to make radio measurements of the flux densities of these components in several different ways. The first and most obvious involves making radio maps of similar resolution to the BIMA data (as in Fig. 3), and then integrating over fixed regions. In addition to the hotspots, we have measured a flux density from the E lobe of the source, weakly detected in the BIMA observations. The regions used are shown on Fig. 4. Secondly, we can also fit a Gaussian and background to identical regions centred on the hotspots. Both integration and Gaussian fitting will tend to produce values of radio flux density which include a contribution from the extended emission around the hotspot.

Thirdly, we can compare the 85-GHz flux densities to the radio flux densities of the most compact radio components only. This was the approach taken by Meisenheimer et al. (1989) in their study of the W hotspot. The assumption here is that the most compact component will be flat-spectrum and will dominate increasingly at higher frequencies. For this purpose we take flux density values for the compact hotspots at 4.9 and 15 GHz from H94, and measure values at 8.5 GHz by fitting a Gaussian and baseline to the full-resolution map of Fig. 1. This technique can be applied only to the W and NE hotspots. There is some ambiguity in applying it to the NE hotspot, which at high resolution has double structure (a bright component and a weak component separated by 0.3 arcsec), as shown in Fig. 1. For consistency with the values of H94 we tabulate the flux density of the brighter component only; for reference, the integrated 8.5-GHz flux density of the weaker component is about 23 mJy.

Results from these three approaches are tabulated in Table 1. The compact components in the NE and W hotspots are well fitted by power laws with spectral indices of $`0.72\pm 0.04`$ and $`0.63\pm 0.02`$ respectively – these values agree well with the two-point spectral indices of H94. There is no evidence for a spectral break at cm wavelengths as there was in the hotspots of 3C 123 (Looney & Hardcastle 2000).

In most cases flux densities measured by integrating regions of the source at 3-arcsec resolution are similar to fluxes derived from Gaussian fitting, as expected. In the remaining analysis we will only use the fluxes measured by integration.

### 4.1 The eastern double hotspot

It is immediately obvious from Table 1 that the overall spectral indices of the NE and SE hotspot components are quite similar; the spectral index over the full frequency range, $`\alpha _{1.4}^{85}`$, is $`0.83\pm 0.03`$ for the NE hotspot and $`0.86\pm 0.02`$ for the SE hotspot. The two spectra are plotted in Fig. 5. There is therefore no immediate evidence for differential spectral ageing in the more extended SE component. Could this just be due to contamination of the NE hotspot flux densities by steep-spectrum extended material included in the integration region? To test this, we subtracted the best-fit power-law spectrum of the compact components from the integrated flux densities, assuming that both compact components of the NE hotspot have the same power-law spectral index of 0.72. The resulting ‘extended’ flux densities for the NE hotspot are also tabulated in Table 1, and it can be seen that almost all the 85-GHz flux of the hotspot can be accounted for as an extrapolation of the spectrum of the compact components. The spectrum of the extended component of the NE hotspot is plotted in Fig. 5. The integrated spectrum of the NE hotspot can therefore consistently be modelled as a sum of compact, flat-spectrum ($`\alpha =0.72`$) and extended, steep-spectrum ($`\alpha 1`$) components, with the expended component possibly cutting off before 85 GHz. But it is clear that the spectrum of the SE component is consistently flatter than the spectrum of the extended component of the NE hotspot, except at low energies where they are comparable.

In detail, the spectrum of the SE hotspot is itself somewhat peculiar; the spectral index in the range 5–15 GHz is around $`0.93\pm 0.01`$, consistent with H94’s estimate of $`0.91`$ with a scatter of $`0.08`$, but it then seems to flatten slightly, with $`\alpha _{15}^{85}=0.89\pm 0.03`$. The inconsistency is only marginal, but it may indicate that a multi-component model is also necessary to explain the SE hotspot’s spectrum; perhaps the SE hotspot is also a superposition of flat- and steep-spectrum components. If so, it is puzzling that H94 find a relatively uniform spectral index across the hotspot, a result confirmed by our own maps of $`\alpha _{1.4}^{8.5}`$ (Fig. 4); there is no obvious site for a flat-spectrum region.

What are the implications of these results for models of the double hotspot? We begin by modelling the spectrum of the compact component of the NE hotspot. Using an angular size derived from our fits to the 8.5-GHz map (which are consistent with the fits of H94), we can model the hotspot as a sphere of radius 380 pc. Assuming that the electron distribution follows a power law with a minimum energy of $`5\times 10^7`$ eV, that the electron power-law index is 2.44 (equivalent to the observed spectral index of 0.72), and that no protons are present in the hotspot, we obtain an equipartition222We feel justified in using equipartition arguments to estimate the field strengths in hotspots, given the detections of X-ray emission at levels consistent with inverse-Compton emission at equipartition in the hotspots of Cygnus A (Harris, Carilli & Perley 1994), 3C 295 (Harris et al. 2000) and 3C 123 (Hardcastle et al. in prep.). field strength fit of 36 nT. Using this, we can constrain the high-energy cutoff of the electron spectrum. The maximum value of the high-energy cutoff consistent with the upper limits on optical flux density derived by H94 is $`1.6\times 10^{11}`$ eV; the minimum possible high-energy cutoff if this component is to account for most of the 85-GHz flux density from the NE hotspot is around $`9\times 10^9`$ eV ($`3\sigma `$ lower limit).

The secondary hotspot’s spectrum is not very consistent with an aged-synchrotron model, because of the flattening of the spectrum between 15 and 85 GHz. It is poorly fit even with a power-law model (best-fit $`\chi ^2=14`$ with 3 degrees of freedom for $`\alpha =0.84\pm 0.01`$, where the errors include a contribution from the uncertainty in the flux scales). If we take the slope of this fit as an indication of the ‘injection index’, the spectral index corresponding to the low-energy electron energy distribution power-law index, and we assume standard Jaffe & Perola (1973, hereafter JP) ageing in a 6-nT field (the equipartition value), then the 85-GHz data point allows us to rule out at the 99 per cent confidence level ($`\mathrm{\Delta }\chi ^2=6.6`$) ages for the SE hotspot greater than $`4\times 10^4`$ years; in other words, particle acceleration in the hotspot is either still going on (in which case a JP spectrum is not appropriate) or it ceased less than $`4\times 10^4`$ years ago. This gives a disconnection timescale of less than one per cent of the total age of the source inferred by spectral ageing methods from the spectrum of the oldest material (e.g. Stephens 1987) but, as we pointed out in our discussion of 3C 123 (Looney & Hardcastle 2000) numerical simulations show that hotspots are transient features on this sort of timescale. We note that if particle acceleration is going on in the hotspot, then this calculation does not tell us anything useful about the hotspot’s age. In this case we would expect a continuous injection (CI) spectrum rather than a JP spectrum (e.g. Pacholczyk 1970), and the position of the break in CI spectra does not bear the same simple relation to the ‘age’ of the electron population as the break in JP spectra.

In the next section we use this spectral information to discuss models for the formation of the double hotspot.

### 4.2 Models for the double hotspot

We can rule out what might be called a trivial dentist’s-drill model in which the jet, having made a single NE-type hotspot at the SE hotspot position, moves to the NE position and makes another, leaving the SE hotspot to expand. In this respect 3C 20 is an extreme example of the arguments put forward by e.g. Valtaoja (1984). The SE hotspot is not only slightly more luminous, but also very much larger than the NE hotspot; for the purposes of modelling, we treat it as a hemisphere of radius 2 arcsec (8 kpc). The ‘trivial dentist’s-drill model’ is then not supported by the energetics of the situation; the (equipartition) total energy in the SE hotspot is $`10^2`$ times greater than that in the NE hotspot (a point which we return to below). But it is even more strongly ruled out by expansion losses, without the need for the assumption of equipartition. To get from a component comparable to the NE hotspots to the SE hotspot by linear expansion, the expansion factor would need to be $`20`$. Adiabatic losses coupled with magnetic field conservation (e.g. Longair, Ryle & Scheuer 1973) would mean that the radio flux density of the progenitor of the SE hotspot before expansion would have been absurdly high, a factor $`10^6`$ higher than its present value; but they would also mean that the high-energy cutoff in the hotspot would now be $`8\times 10^9`$ eV, based on the upper limit on the high-energy cutoff in the NE hotspot, which, given an equipartition field strength of $`6`$ nT estimated from the observations, is far too low to produce the observed 85-GHz emission from the SE hotspot. The secondary component in 3C 20 cannot have evolved by adiabatic expansion from a single primary with properties even slightly resembling those of the NE hotspot, given a simple geometry. The problem is only slightly alleviated if we consider less simple geometries for the emitting material in the southern hotspot.

Instead, we can consider a more sophisticated version of the dentist’s drill model, in which the jet termination point moves around in the region of the SE hotspot before switching to that of the current primary. In this case, the energy available to the SE hotspot is potentially that of a large number of components more similar to the NE hotspot, which reduces the problem of energetics. However, the adiabatic expansion argument is harder to evade in this way; the multiple versions of the primary hotspot left in the volume of the SE hotspot will all expand adiabatically once the jet moves on to the next one, and so each in turn will become fainter by a large factor. Clearly this model requires just on energetic grounds that the jet has dwelt in the region of the SE hotspot at least a hundred times as long as it has been powering the NE hotspot, and probably much longer once adiabatic losses are taken into account. For this model to be viable we therefore have to be seeing the source at a very special time in its history, and while this is not impossible it is not an explanation that we wish to adopt while others exist.

The other possible models all require that there is, or was until recently, long-term energy input into the SE hotspot from a source other than a ‘primary-type’ jet. The energy input may either be simple advection of pre-accelerated particles, or it may predominantly be in the form of kinetic energy of a beam, with particle (re)acceleration taking place in the SE hotspot. We comment on several possible models in turn.

* Disconnected-jet models. In these models, described by Cox et al., the remains of a disconnected jet continue to power the secondary hotspot after a new primary has started to form. This would help to explain the low age of the material in the secondary hotspot. We see two problems with these models in the case of 3C 20. One is that the relaxed appearance of the hotspot suggests that there is not currently well-collimated energy input – a disconnected jet’s impact point should presumably look similar to a primary hotspot’s, since the material of the disconnected jet presently entering the hotspot does not ‘know’ that it has been cut off further upstream. The other is that this model does not help to account for the factor of 100 difference between the energy content of the SE and NE hotspots; like the dentist’s drill model, it requires that we are seeing the source at a special time.
* Two-jet models. We cannot rule out the possibility that there are two jets, or a split jet, with a powerful but poorly collimated jet powering the secondary while a weak, well-collimated jet powers the primary. The only evidence against this is that the ‘secondary jet’ this would require has never yet been observed in any radio source.
* Outflow from the primary. Laing (1982) was the first to point out that the SE hotspot is brightest in the region opposite the compact NE hotspot. That suggests that the secondary may be powered by direct outflow from the primary, as in the splatter-spot or redirected-outflow models. The light travel time from the primary to the secondary hotspot is $`10^5`$ years (the uncertainty is due to the unknown projection angle). So the age constraints derived above rule out simple advection of particles from the primary, unless they are transported at very high speeds ($`\gamma _{\mathrm{bulk}}2`$) or in a low-loss channel. Instead, it seems most likely that energy is transported to the SE hotspot largely as bulk kinetic energy, in quite poorly collimated outflow from the primary, and that particle acceleration is ongoing there. (The flattening of the secondary hotspot spectrum between 15 and 85 GHz, and the fact that the secondary hotspot’s spectrum is flatter than that of the extended component of the primary, may provide some evidence for particle acceleration.) In this model, the energy content of each hotspot is a function not just of their history, but of the rate at which energy is transported out of the primary and of the efficiency of the particle acceleration process (there is no reason in principle why the primary cannot be significantly less efficient than the secondary in turning the bulk kinetic energy provided by the jet into internal energy of relativistic particles). The large difference between the energies of the two hotspots therefore does not need to be explained in terms of a special observing time in this picture.

We conclude that the model with the fewest objections against it is the outflow model, and therefore join Cox et al. (1991) in favouring that model for this source.

### 4.3 The western hotspot

The integrated spectrum of the W hotspot (Fig. 6) is not well modelled as a simple power-law. Since the region of the BIMA detection contains components with a wide range of spectral indices (Fig. 4), this is probably not surprising. Table 1 shows that the most compact component contributes less than 20 per cent of the total flux density of the integration region even at 15 GHz. The northern ‘tail’ of the compact hotspot (Fig. 1) contributes slightly less, about 150 mJy at 8.5 GHz. If we assume the tail and hotspot have the same power-law spectral index of $`0.63\pm 0.02`$ (H94 suggest a steeper spectrum for the northern tail, but this may in part be an effect of a steep-spectrum background) then we can subtract off an extrapolation of the compact components to obtain the estimate of the ‘extended’ spectrum tabulated in Table 1 and plotted in Fig. 6. Like the extended spectrum of the NE hotspot, this is not well fitted either by a power-law model or by an aged-synchrotron model. But the spectral shapes of the W and NE extended components are similar, which suggests that the same sort of physical processes (perhaps related to escape from the compact hotspots) are at work.

### 4.4 The lobes

Fig. 4 shows that the lobes of 3C 20 have steep spectra, which is conventionally attributed to spectral ageing. Some authors have questioned whether the spectra of lobes can really be described in this way. Katz-Stone, Rudnick & Anderson (1993) show that colour-colour diagrams (where the spectral indices for one pair of frequencies are plotted against the spectral indices for another) are useful indicators of whether sources follow any of the standard ageing models, which correspond to well-defined lines in the colour-colour plane; they argue that colour-colour diagrams for Cygnus A are not consistent with any of the standard ageing models. Rudnick (1999) gives a summary of the problems of spectral ageing.

In Fig. 7 we plot a colour-colour diagram for 3C 20, using the 1.4-, 4.9- and 8.5-GHz VLA data, together with lines showing the tracks of JP aged synchrotron spectra with injection indices of 0.5 and 0.75. The results are quite similar to those of Katz-Stone et al. The bulk of the source material can be adequately modelled as a JP aged synchrotron spectrum (with an injection index $`0.75`$) but at flat spectral indices the source deviates from the expected spectrum; we do not see points which are described by an un-aged power law with $`\alpha =0.75`$. In 3C 20, unlike Cygnus A, some source material is well described as a simple power law, but this has $`\alpha 0.6`$. Breaking down the points in Fig. 7 by source region, we find that the points with $`\alpha _{1.4}^{4.9}<0.8`$ and $`\alpha _{4.9}^{8.5}<0.9`$ are all in the hotspot regions shown in Fig. 4, with the $`\alpha 0.6`$ points being contributed by the W hotspot. The points with $`\alpha _{1.4}^{4.9}>0.8`$ and $`\alpha _{4.9}^{8.5}>0.9`$, which can be modelled with a JP spectrum, are essentially all in the lobes of the source. The observation that there are no points with $`\alpha _{1.4}^{4.9}\alpha _{4.9}^{8.5}0.75`$ thus corresponds to the observation that spectral ages in the lobes of other FRII sources often do not extrapolate back to zero age as distance from the hotspots tends to zero (Alexander 1987, Alexander & Leahy 1987). Stephens (1987) found a similar effect in a low-resolution spectral age study of 3C 20.

The outstanding questions raised by Fig. 7 are:

1. Is the apparently adequate description of the lobes with JP aged synchrotron spectra indicating the real underlying physics, or is it a coincidence? If the lobes really are simply aged material,
2. why is the injection index of the lobes apparently different from the (flatter) injection indices in the compact hotspots, and
3. why do we see no source material which has a simple power law spectrum with spectral index equal to the injection index; equivalently, why is there no zero-age material with the injection index which appears to characterize the lobes?

Since the hotspots are transient features, question (ii) can be avoided, if necessary, by noting that the injection index of the present-day hotspots need not be the same as the injection indices of the hotspots present at the time that the lobe material was accelerated. Alexander (1987) suggests that the answer to question (iii) may be that the effects of losses in the comparatively high field of the hotspot are exaggerated by adiabatic expansion, which for linear expansion by a factor $`\mathrm{\Delta }`$ moves the break frequency down by a factor $`\mathrm{\Delta }^4`$ if magnetic flux is conserved. For $`\mathrm{\Delta }2`$ as suggested by Alexander in the case of 3C 234, this would imply that we should start to see ageing in the mm region of the hotspot spectra, for which there is no convincing evidence in the case of 3C 20; however, because of the low resolution of our BIMA observations, we cannot really rule out the presence of steeper-spectrum regions around the flat-spectrum hotspots.

An answer to question (i) is beyond the scope of this paper, but we note that the observation of 85-GHz emission from the lobe region is very consistent with the predictions of a JP aged model (Fig. 8); if the similarity to a JP spectrum is a coincidence, it is one which extends over a wide range of frequencies.

## 5 Conclusions

BIMA observations of the double hotspot of 3C 20 resolve the two components of the hotspot and, together with archival VLA data, show that their radio spectra are very similar over the full range studied; 3C 20 in this respect is similar to other double hotspot sources like 3C 405 and 3C 123 previously studied with BIMA. The detailed spectrum of the secondary hotspot in 3C 20 provides some evidence that particle acceleration is continuing, and this, together with the relaxed and edge-brightened appearance of the secondary in high-resolution radio maps, causes us to conclude that the most likely model for this source is one in which the secondary is generated by direct outflow from the primary hotspot; this outflow can be relatively slow and poorly collimated compared to the jet.

Both the eastern and western compact (primary) hotspots are bright sources at 85 GHz, and it seems likely that there is some emission at this frequency from the regions around them as well. Our radio observations agree with H94’s statement that the radio spectrum of the eastern primary is slightly but significantly steeper than that of the W primary, but the 85-GHz observations are consistent with the primary hotspots being otherwise similar at high frequencies; in fact, our estimated upper limit on the high-energy cutoff in the compact component NE hotspot is exactly consistent with the high-energy cutoff required to produce the observed optical emission from the compact component of the W hotspot, if the flatter injection index and slightly higher equipartition field strength in the W hotspot are taken into account. Further observations at frequencies intermediate between 85 GHz and optical will be required to establish whether the primary hotspots, so similar in size and energy content, have different high-energy cutoffs in their electron spectra, and to establish a difference between the NE and SE hotspots.

Our VLA and BIMA data show that 3C 20’s lobes are quite well described by a simple JP ageing model, but that, as with other sources, it is difficult to understand either why the injection indices in the lobes and hotspots are different or why no material with a pure injection-index spectrum is seen.

## ACKNOWLEDGEMENTS

We thank an anonymous referee for comments which enabled us to improve the paper significantly. We are grateful to Rick Perley for allowing us to use the archival 5- and 15-GHz VLA data discussed in the text, and the VLA data analysts for help in extracting it from the archive. The National Radio Astronomy Observatory VLA is a facility of the U.S. National Science Foundation operated under cooperative agreement by Associated Universities, Inc.


text:
# The stress field in granular systems: Loop forces and potential formulation

## Abstract

The transmission of stress through a marginally stable granular pile in two dimensions is exactly formulated in terms of a vector field of loop forces, and thence in terms of a single scalar potential.

The loop force formulation leads to a local constitutive equation coupling the stress tensor to fluctuations in the local geometry. For a disordered pile of rough grains (even with simple orientational order) this means the stress tensor components are coupled in a frustrated manner, analogous to a spin glass.

If the local geometry of the pile of rough grains has long range staggered order, frustration is avoided and a simple linear theory follows. Known exact lattice solutions for rough grains fall into this class.

We show that a pile of smooth grains (lacking friction) can always be mapped into a pile of unfrustrated rough grains. Thus it appears that the problems of rough and smooth grains may be fundamentally distinct.

It has long been recognized in engineering practice that granular materials which lack cohesion cannot be regarded entirely as solids, nor (until ’fluidised’ or yielding) can they be regarded entirely as fluids . In this Letter we explore the scenario that there exists a marginal state of granular matter in between solid and liquid. Most particularly, we are concerned with how that intermediate state transmits stress.

The marginal state is readily characterised in terms of the mean coordination number $`\overline{z}`$ for mechanical equilibrium. If external loads, such as gravity, are exerted on a pile of $`N`$ grains  then the $`N\overline{z}/2`$ intergranular forces must be able to adjust so as to achieve $`Nd(d+1)/2`$ constraints balancing the force and torque on each grain in $`d`$ dimensions. For perfectly rough grains, where every contact supports friction, this gives a critical coordination number $`z_c=d+1`$, whilst for ideally smooth asymmetric grains with no friction $`z_c=d(d+1)`$ ,,.

For mean coordination numbers $`\overline{z}<z_c`$ the pile cannot be stable and under a general loading it must rearrange or consolidate. By contrast for $`\overline{z}>z_c`$ the intergranular forces are underdetermined by the conditions of force balance alone and the deformation of individual grains together with their local constitutive equations becomes relevant. In the first case we have granular motion and the pile is fluid or yielding, whereas in the latter case intragranular deformation is crucial and we have a solid of constitution influenced by that of the individual grains. It is interesting to note that for ideally smooth grains, $`z_c`$ matches the maximum mean coordination attainable by sequential packing (and in two dimensions it is the absolute maximum possible ), so a (random) pile of smooth grains is probably never solid.

The marginal case $`\overline{z}=z_c`$ is special in that the grains do not need to move and the intergranular forces, and hence distribution of stress, are determined by conditions of force balance alone. This has been identified as a paradigm problem of theoretical granular mechanics , , . The macroscopic analogues of balancing force and torque are that the stress tensor $`\sigma `$ obey

$$\stackrel{}{}\widehat{\sigma }=\stackrel{}{F},𝒜(\widehat{\sigma })=0,$$
(1)
where $`𝒜(\widehat{\sigma })=\frac{1}{2}(\widehat{\sigma }\widehat{\sigma }^T)`$ is the antisymmetric part of $`\widehat{\sigma }`$ and $`\stackrel{}{F}`$ is an external force field. We shall focus here on a region transmitting stress but otherwise unloaded. This system of continuum equations is incomplete without $`d(d1)/2`$ further constitutive relations. The central mystery we are seeking to unravel is the nature of these relations (one relation in two dimensions) when they also arise from force balance at the granular level.

The key idea of this Letter, where we limit ourselves for simplicity to the two dimensional case, is to focus attention on local loops around contact points and corresponding loop forces. Each local void, which we will label by an index $`l`$, is enclosed by a loop of grains and contact points; each contact point separates two voids and so is party to two loops. Around each loop we take a loop force $`\stackrel{}{f}_l`$ to circulate in the anticlockwise direction, that is $`+\stackrel{}{f}_l`$ is contributed to each intergranular force around the loop with (notionally) positive sign in the anticlockwise direction, and vice versa for the reaction forces. The resulting total force across a contact point is then a difference between two contributing loop forces.

The loop forces have two key features, the first being that they parameterise the intergranular forces in a manner that automatically satisfies balance of force at each contact point and on each grain. In terms of them we can write the stress tensor on each grain as

$$\widehat{S}_g=\underset{l}{}\stackrel{}{r}_{lg}\stackrel{}{f}_l$$
(2)
where $`\stackrel{}{r}_{lg}`$ is the vector connecting (anticlockwise) the two contact points shared by loop $`l`$ and grain $`g`$ and $`A_g`$ is the grain area defined in fig. 1 . The vectors $`_g\stackrel{}{r}_{lg}=0`$ form an anticlockwise loop round the contact points of loop $`l`$ and likewise $`_l\stackrel{}{r}_{lg}=0`$ forms a clockwise loop round the contact points of grain $`g`$ and play a central geometrical role in our discussion.

The second key feature of the loop forces is that, being defined on the loops which are half as numerous as the grains, they are a comparatively coarse-grained quantity. The hidden nature of the constitutive equation can now readily be appreciated, because the only remaining constraints that the loop forces have to obey is balance of torque around each grain, that is

$$\underset{l}{}\stackrel{}{r}_{lg}\times \stackrel{}{f}_l=0$$
(3)
which is one equation per grain and hence two equations per loop.

In essence the two torque equations per loop give us two macroscopic conditions, $`𝒜(\widehat{\sigma })=0`$ plus a constitutive relation. To achieve this separation explicitly we postulate a smooth interpolation of the loop forces to a function of continuous position $`\stackrel{}{f}(\stackrel{}{r})`$, with each loop and grain having a nominal centre $`\stackrel{}{R}_l`$ and $`\stackrel{}{R}_g`$, respectively. This also introduces a natural continuum already at the grain scale without going into elaborate coarse-graining schemes. Then the force moment on grain $`g`$ is, to a first order Taylor approximation,

$$\widehat{S}_g=\widehat{C}_g\stackrel{}{}\stackrel{}{f}$$
(4)
where

$$\widehat{C}_g=\underset{l}{}\stackrel{}{r}_{lg}\stackrel{}{R}_{lg}\text{and}\text{ }\stackrel{}{R}_{lg}=\stackrel{}{R}_l\stackrel{}{R}_g.$$
(5)
The geometric tensors $`\widehat{C}_g`$ characterise the grain local geometry and their symmetric and antisymmetric parts have quite distinctive properties. It is readily shown that

$$\widehat{C}_g=A_g\widehat{R}+\widehat{P}_g$$
(6)
where $`\begin{array}{c}\widehat{R}\end{array}=\left(\genfrac{}{}{0pt}{}{\mathrm{0\hspace{0.17em}1}}{\mathrm{1\hspace{0.17em}0}}\right)`$ is the unit antisymetric tensor (also corresponding to $`\frac{\pi }{2}`$-rotation in the plane) and $`A_g`$ is the area associated with grain $`g`$ shaded in figure 1. The areas $`A_g`$ have the convenient property of tessalating the plane. The symmetric part, $`\widehat{P}_g`$, reflects fluctuation in the departure of the local geometry from isotropy, being given by

$$\widehat{P}_g=\underset{l}{}\frac{1}{2}(\stackrel{}{s}_{lg}\stackrel{}{s}_{lg}\stackrel{}{t}_{lg}\stackrel{}{t}_{lg})$$
(7)
where the vectors $`\stackrel{}{s}_{lg}`$ and $`\stackrel{}{t}_{lg}`$ are shown on figure 1. It is a crucial point that neighbouring grains share equal and opposite contributions $`\stackrel{}{s}\stackrel{}{s}`$ and $`\stackrel{}{t}\stackrel{}{t}`$ to their $`\widehat{P}`$ tensors, and from this it follows that the average $`\widehat{P}_G`$ of $`\widehat{P}`$ over all the grains in any convex region $`G`$ goes to zero at least as fast as $`N_G^{1/2}`$, where $`N_G`$ is the number of grains in $`G`$. Note that the averaging of the $`P_g`$ does not depend on any assumptions of isotropy. For periodic lattices $`\widehat{P}`$ vanishes over one unit cell.

To underpin the significance of $`C_g`$ as the relevant geometrical order parameter, consider a region of area $`A`$ whose boundary is defined by the contact points of the outermost grains. Averaging over $`\widehat{S}_g`$ in the region using eq. (4) gives

$$\widehat{\sigma }=\frac{1}{A}\underset{g}{}\widehat{C}_g\stackrel{}{}\stackrel{}{f}=\widehat{C}_g\stackrel{}{}\stackrel{}{f}.$$
(8)
Carrying out the same summation using eq. (2) we observe that the sum over loops inside the region cancels out and the only contribution comes from the boundary vectors, $`\stackrel{}{r}_b`$, between boundary contact points,

$$\widehat{\sigma }=\frac{1}{A}\underset{\mathrm{boundary}}{}\stackrel{}{r}_b\stackrel{}{f}_b,$$
(9)
where $`\stackrel{}{f}_b`$ is the external loading on the region. This contour sum can be converted, using Stokes theorem, into

$$\widehat{\sigma }=\frac{1}{A}𝑑s\widehat{R}\stackrel{}{}\stackrel{}{f}=\widehat{R}\stackrel{}{}\stackrel{}{f}.$$
(10)
Since $`\widehat{P}_g=0`$ we identify $`\widehat{C}_g=A\widehat{R}`$ and, on comparing to eq. (8), we obtain

$$\widehat{C}_g\stackrel{}{}\stackrel{}{f}=\widehat{C}_g\stackrel{}{}\stackrel{}{f}.$$
(11)
This result is particularly reassuring because it yields the effective response of a macroscopic granular region to a force field in terms of an effective macroscopic characteristic property $`\widehat{C}_g`$. It is equivalent to results in contexts such as disordered dielectrics and continuum elasticity that are derived by combining a field equation with either a constitutive relation or an energy functional.

We can now construct the constitutive equation in terms of the continuum stress tensor, which is given from eq. (10) as $`\widehat{\sigma }=𝒟\stackrel{}{f}`$ where $`𝒟`$ is exactly the curl operator in two dimensions. The condition that the stress tensor on each grain be symmetric then leads, over and above the condition that the continuous $`\widehat{\sigma }`$ be symmetric, to the requirement

$$\mathrm{Tr}(\widehat{Q}\widehat{\sigma })=0,$$
(12)
where $`\widehat{Q}=\widehat{R}\widehat{P}\widehat{R}^T`$ is a rotated verion of the symmetric part of the local $`\widehat{C}_g`$ tensor. Relation (12) is the continuum constitutive equation that provides the missing link between the stress and the local geometry. It has the striking feature that the coefficients $`\widehat{Q}`$ are spatially fluctuating quantities that locally add to zero, so that any simple attempt to identify a non-vanishing mean field value $`\widehat{Q}`$ leads back to Eq. (1) and therefore yields no new constitutive information. In this sense the problem is analogous to spin glasses, in our case it being $`\widehat{\sigma }`$ that is subject to spatial couplings of random sign .

Constitutive equations coupling stress components in a local linear relation, represented here by (12), have been discussed before under the presumption that a non-vanishing macroscopic version of $`\widehat{Q}`$ exists. Examples are the hypotheses of ”fixed principal axes”  (corresponding to $`\widehat{Q}`$ being traceless) and ”oriented stress linearity”  that were introduced to describe sand piles . Having shown that the simple grain average of $`\widehat{Q}`$ vanishes, we are therefore particularly motivated to consider any case where a non-zero effective $`\widehat{Q}`$ can be identified.

Drawing inspiration from spin problems, we focus on the analogy of an antiferromagnet in which the local $`\widehat{P}`$ tensors have alternation of sign in an unfrustrated pattern. The key idea is that for systems where every loop has an even number of edges, we can label all grains \+ or - such that each grain is surrounded only by opposite sign neighbours in an unfrustrated way. Then we can distinguish locally averaged values $`\widehat{P}_+`$ and $`\widehat{P}_{}`$ restricted to the respective $`+`$ and $``$ grains, where $`\widehat{P}_++\widehat{P}_{}=0`$. If the pile has staggered geometrical order such that these are non-zero, then this leads to the identification

$$\widehat{Q}_{\mathrm{eff}}=\widehat{R}\widehat{P}_+\widehat{R}^T$$
(13)
The simplest such case where we can see this work is a periodic lattice (fig. 2) with a generally anisotropic unit cell comprising of two grains, which we duly label $`+`$ and $``$. Using eq. (4) for each grain separately, we write down the mean and deviatory cell stresses, for which we find

$$\widehat{\sigma }=\left(\widehat{S}_++\widehat{S}_{}\right)=\left(A\widehat{R}+\widehat{P}_++\widehat{P}_{}\right)\stackrel{}{}\stackrel{}{f}$$
(14)
and

$$\widehat{\delta }=\left(\widehat{S}_+\widehat{S}_{}\right)=\left(\delta A\widehat{R}+\widehat{P}_+\widehat{P}_{}\right)\stackrel{}{}\stackrel{}{f},$$
(15)
where $`A`$ is the unit cell area and $`\delta A=A_+A_{}`$ is the difference between the areas associated with the two grains. Using the periodic properties of the lattice one readily obtains

$$\widehat{P}_+=\widehat{P}_{}=\frac{1}{2}\underset{k}{}\stackrel{}{s}_k\stackrel{}{s}_k\stackrel{}{t}_k\stackrel{}{t}_k$$
(16)
where the vectors $`\stackrel{}{s}_k`$ and $`\stackrel{}{t}_k`$ are labelled on figure 2. Then the condition that both $`\widehat{\delta }`$ as well as $`\widehat{\sigma }`$ be symmetric then leads to the constitutive equation (12) with $`\widehat{Q}`$ replaced according to (13).

The idea of a staggered order parameter can be further exploited to formulate a mapping from perfectly smooth to perfectly rough systems. Given a pile of perfectly smooth grains, we note that its behaviour is identical to the same pile of rough grains provided at each contact point we insert a vanishingly small perfectly rough ball bearing. It is readily checked that starting from $`\overline{z}=6`$ for the smooth grains leads to $`\overline{z}=3`$ for the resulting set of rough grains (3/4 of which are the new bearings each having $`z=2`$), so we duly arrive at a marginally connected set of rough grains. Furthermore, we automatically have an alternating distinction between original grains ($`+`$) and bearings ($``$). Explicit evaluation shows that for each small grain, $`\widehat{P}_{}ϵ`$ as bearing radius $`ϵ0`$, so in this case the appropriate geometrical order parameter is $`p=Limit_{ϵ0}P_+/ϵ`$. This order parameter, and hence the corresponding $`Q_{\mathrm{eff}}`$, can have a non-zero average for an orientationally ordered assembly of smooth grains, as computed explicitly for a periodic array. This contrast to the case of general rough grain piles, which explicitly require staggered orientational order for a non-zero $`Q_{\mathrm{eff}}`$, suggests that piles of perfectly smooth grains may be inherently distinct from piles of rough ones.

Finally we introduce an alternative formulation of the problem where the loop forces are defined in terms of a potential field. Since each loop is associated with a two-dimensional force the number of degrees of freedom that the potentials should provide is twice the number of loops, which (up to relative correction of $`1/\sqrt{N}`$ due to the boundary) is exactly the number of grains in the system. Thus, we require that the potential field provide one degree of freedom per grain leading to scalar grain potentials, $`\psi _g`$. A natural way to precisely define the potentials is to choose that the forces are derived from them via

$$\stackrel{}{f}_l=\frac{1}{A_l}\underset{g}{}\stackrel{}{r}_{lg}\psi _g,$$
(17)
and the potentials are then governed by eq. (3).

We now assume that the $`\psi `$-field can be continued throughout the system just as we did previously for the force field. Then expanding around the loop centres, we obtain

$$\stackrel{}{f}_l=\frac{1}{A_l}\widehat{C}_l\stackrel{}{}\psi .$$
(18)
The geometrical tensor $`\widehat{C}_l=_g\stackrel{}{r}_{lg}\stackrel{}{R}_{lg}`$ is the loop analogue of $`\widehat{C}_g`$ and its antisymmetric part is likewise associated with an area $`A_l`$, where these areas tessalate the plane. The symmetric part of $`\widehat{C}_l`$ is sensitive to the choice of grain centre positions: as there are two grains per loop these constitute four degrees of freedom per $`\widehat{C}_l`$ tensor and the expectation is that these tensors can thereby be chosen to be purely antisymmetric. Using $`\widehat{C}_l=A_l\widehat{R}`$ simplifies the expression of the local stress in terms of the potentials to

$$\widehat{S}_g=\widehat{C}_g\widehat{R}^T𝒟_g𝒟_l\psi $$
(19)
and the continuum stress becomes

$$\widehat{\sigma }=𝒟_g𝒟_l\psi \mathrm{𝐜𝐮𝐫𝐥}\mathrm{𝐜𝐮𝐫𝐥}\psi .$$
(20)
It follows that the continuous version of $`\psi `$ is exactly the Airy stress function .

To conclude, we have formulated a theory of stress transmission in granular matter at the marginal state between solid and fluid in terms of loop forces, which in turn can be defined in terms of a scalar potential. Granular systems can possess a predominant antiferromagnetic order but with local frustration, much like in antiferromagnetic Ising systems. An antiferromagnetic order yields a global consitutive relation leading to a macroscopic theory while frustration requires a more sophisticated approach. The similarity between random granular piles and spin glasses suggests a possible use of methods from the latter to this problem. We have mapped smooth piles to antiferromagnetically ordered perfectly rough ones obtaining for them an exact constitutive relation. This result points to an inherent difference between rough (frustrated) and smooth (ordered) systems and may thus undermine the idea of gaining insight into behaviour of piles of finite friction grains from understanding the perfectly smooth and perfectly rough limits.

An extension of this work to three dimensions is currently under way.

Acknowledgements

We are grateful to Prof Sir S. F. Edwards and Dr D. Grinev for many illuminating discussions. This research was supported by EPSRC grant GR/L59757.