text:
# Multispeckle diffusing-wave spectroscopy: a tool to study slow relaxation and time-dependent dynamics

## I Introduction

Dynamic light scattering has proven to be a powerful tool for investigating the dynamics of colloidal suspensions, emulsions, foams, and other complex fluids. Over the past decade, dynamic light scattering (DLS) has been increasingly applied to glassy systems with slow, and often non-ergodic, dynamics. Slow dynamics pose experimental difficulties for DLS, particularly with respect to the stability of the lasers, and sometimes for the light detection system. Nonergodic samples pose even more serious experimental difficulties, which are exacerbated by the presence of slow dynamics. One way of addressing the problem of long data acquisition times is to use an array detector. More importantly, an array detector can also address the problem of non-ergodic samples if the optics are designed so that the different elements in the array detector constitute essentially statistically independent samplings of the system under study.

In this article, we describe a method which employs a charge-coupled device (CCD) array detector for performing dynamic light scattering on systems that exhibit strong multiple light scattering. While the basic method we describe is applicable to single as well as multiple light scattering, our focus here is on multiple dynamic light scattering, also known as diffusing-wave spectroscopy (DWS) maret87 ; pine88 ; weitz93 , where the methods we describe are most effective.

## II Background

In conventional DWS (or DLS), the basic idea is to illuminate a sample with a coherent beam of light from a laser and to measure the temporal fluctuations in the resulting speckle pattern of scattered light. In a typical experiment, the optics for collecting the scattered light is arranged in such a way that the intensity of a single speckle is measured. By focusing on a single speckle, the amplitude of the fluctuating intensity is optimized to obtain a strongly fluctuating signal. The temporal duration of a typical fluctuation, i.e. the lifetime of a speckle, is determined by the time it takes a particle to move such that the phase of the scattered light is changed by approximately unity (or $`\pi `$). The dynamics of the scattered light, and hence the dynamics of the scatterers, can be more quantitatively characterized by the temporal autocorrelation function of the scattered intensity. In many cases, the autocorrelation function can be analyzed so that the time dependent mean square displacement of the scatterers is extracted weitz93 .

There are two important keys to the success of the conventional analysis that we wish to focus on here. First, the process must be stationary, at least for the duration of the experiment. That is, the dynamics of the system cannot change significantly during this time period. For example, the gelation process of a colloidal suspension can be monitored provided that the particle dynamics slow down (due to the increase of the elastic modulus) occurs on time scales longer than the measurement period. In general an accurate measurement of the correlation function requires that the fluctuations in the intensity be measured over many fluctuation periods (e.g. $`10^4`$ speckle lifetimes to achieve $`1`$% accuracy). Thus, in this case, the particle dynamics must not change significantly over a time period of $`10^4`$ speckle lifetimes.

Second, the system must be ergodic. For a light scattering experiment, a system is ergodic if each speckle samples the same distribution of light intensities as the entire ensemble of speckles over the period of time that data is acquired (i.e. over the duration of the experiment). While this condition is satisfied for a large number of systems, there are certain systems with very slow glassy dynamics that are of interest but which, nevertheless, are not ergodic in the sense described above. Such systems pose a problem for light scattering and special steps must be taken in order to assure that meaningful correlation functions, which represent dynamics characteristic of the entire system, are measured pusey89 ; joosten90 ; xue92 .

The problem for light scattering presented by nonergodic samples was first understood in a clear fashion by Pusey and van Megen who, at the same time, presented a clever method of obtaining correctly averaged temporal autocorrelation functions pusey89 . This technique has proven to be very useful but has some drawbacks: it sacrifices some signal-to-noise in order to obtain a useful distribution of light intensities, and the dynamics of the system must not change significantly over the duration of the experiment pusey89 ; joosten90 ; xue92 . Moreover, it suffers from the same limitation as normal DLS experiments in that data must be collected over a time scale that is much longer that the longest time scale of the fluctuations one is trying to measure.

Other schemes have been developed which overcome some of these difficulties, but not without a significant cost. The sample can be translated or rotated, for example, so that an ensemble average is obtained by brute force xue92 . However, the rotation or translation introduces an artificial characteristic time scale into the correlation function, namely the time scale over which the system samples new speckles by virtue of the motion. Care must be taken to assure that this time scale is much longer than any dynamical time scale of the system one is interested in measuring. And once again, the duration of the experiment must be such that averaging occurs over a time scale that is long compared to the time scale introduced by the sample motion (this amounts to assuring that a sufficiently large number of speckles are sampled – e.g. $`10^4`$ – to obtain the desired level of accuracy – e.g. $`1\%`$). Another method, developed recently by Scheffold et al. skipetrov96 ; romer00 ; scheffold01 , uses two cells, one containing the sample with the slow nonergodic dynamics that is being studied, and another containing diffusing particles with ergodic dynamics, which serves to spatially average the nonergodic signal. This technique is only useful for experiments involving strong multiple scattering, but nevertheless provides a simple and efficient means of dealing with nonergodic systems. However, the only time scales probed in the primary cell are those that are shorter than the diffusive time scale set by the secondary cell.

An attractive alternative to the techniques described above for dealing with nonergodic samples is to use a CCD array, conveniently contained in a commercial CCD camera, as a light detector wiltzius93 ; sillescu97 ; tsui98 ; cipelletti99 ; cipelletti00 ; knaebel00 . This allows the intensity from a large number of speckles to be simultaneously and independently measured. A time series of images of scattered light obtained from the CCD camera are stored in a computer and processed such that each pixel in the CCD array serves as an input to its own correlator. Since the fluctuation of each speckle is independent of the other speckles, processing scattered intensity data collected with a CCD camera is like have $`N`$ independent samplings of the data, where $`N`$ is the number of speckles in the field of the CCD array. If $`N`$ is sufficiently large, e.g. $`N10^4`$, then a complete ensemble average is obtained without any need to sample the system over times any longer than the actual dynamical time scales that are under study. Thus, a direct ensemble average can be obtained without any need for time-averaging whatsoever. This approach solves several problems at once. First, it solves the nonergodicity problem because a full ensemble average (of $`10^4`$ samples) is obtained. Second, the process under study need not be a stationary process because no time averaging is necessary, provided a sufficiently large number of speckles are used. Thus, an essentially instantaneous snapshot of the system dynamics can be taken. If those dynamics change in time, they can be followed simply by measuring temporal autocorrelation functions at different starting times.

There are two practical limitations of the CCD-camera-based technique. First, the minimum delay time for which the correlation function can be measured is set by the frame rate of the camera. In the experiments described in this article, we use a frame rate of 500 frames per second. Thus, delay times as small as a few milliseconds are accessible. It is interesting and important to note that real time processing of the images is necessary since, to our knowledge, no commercial buffer memory is able to store them for the entire sequence of frames obtained in the course of a single experiment, which is typically on the order of $`10^5`$ to $`10^7`$ frames. Recent advances in microprocessor speed allows real time computation of the correlation function. Hence the limiting factor for frame rate is the PCI bus velocity for transferring data from the acquisition board to the computer RAM.

The second limitation of the technique is the limited sensitivity of the CCD array used, as compared to that of a photomultiplier. This means that relatively strong scattered light signals are needed. In practice a power of 400 mW gives a satisfactory signal to noise ratio for a frame rate of 500 per seconds.

In this article we first describe the CCD optical set up and calculational scheme for our computer-based correlator. In the next section we explain how a single camera can be used to take up to 500 correlation functions simultaneously, with the reference starting time shifted by different amounts. We then show how the CCD camera technique can be coupled with the double cell technique to extend the time window over which measurements can be made from 10-8 s to 104 s.

## III Experimental Setup

The sample is illuminated using polarized light from an argon-ion laser operating at a wavelength of $`\lambda =514`$ nm. Its maximal output power is 1.6 Watts. The laser beam is expanded to a diameter of approximately 1 cm and is incident on the sample cell, as shown in Fig. 1. Multiply scattered light is collected by a 50-mm Nikon camera lens and focused onto a iris diaphragm. The lens is set up so that a one-to-one image of the scattered light from the output plane of the sample cell is formed on the diaphragm. The imaging of the output plane of the sample onto the pinhole provides a convenient means for collecting data from different parts of the sample. This feature is particularly useful in samples like gels or very viscous solutions where defects or small particles of dust may hinder the measurement; this set-up allows one to avoid defected regions of the sample. A polarizing cube is placed between the lens and the diaphragm to protect the camera from stray light.

The low sensitivity of CCD arrays relative to photomultiplier tubes means that the sample must be strongly illuminated. However, the incident light intensity must not be so high as to heat the sample. Therefore, it is crucial to collect as much scattered light as possible. We use a camera lens of 50 mm with a numerical aperture of 1.7.

The scattered light is detected by a Dalsa CCD camera, model CAD1-128A, which is placed approximately $`d15`$ cm behind the diaphragm. It is an 8-bit camera which we run at 500 frames per second. The images are transferred to a computer running at 500 MHz using a National Instruments data acquisition board, model PCI1422. The combined speeds of the computer and data acquisition boards are sufficiently high to allow data to be analyzed in real time.

Data analysis includes calculation of the ensemble-averaged temporal correlation function of the scattering intensity $`I`$,

$$g_2(t,t_0)\frac{I(t_0)I(t_0+t)}{I(t_0)I(t_0+t)},$$
(1)
and the ensemble-averaged spatial correlation function

$$C(\delta 𝐫)\frac{I(𝐫)I(𝐫+\delta 𝐫)}{I(𝐫)^2}$$
(2)
where $`t_0`$ is the time when the reference intensity is measured, $`tt_0`$) is the delay time, $`𝐫`$ is the coordinate of a pixel on the detector, and $`𝐫+\delta 𝐫`$ is the coordinate of another pixel a distance $`\delta 𝐫`$ away.¡¿ represent the average over the CCD chip. The denominator of the expression on the right hand side of Eq. (1) is the square of the average intensity determined according to the method suggested by Schatzel schatzel91 , in order to obtain a slightly more accurate measurement of $`g_2(t,t_0)`$. The algorithms used to calculate these correlation functions are described in greater detail below.

### III.1 Speckle size adjustment

The size $`s`$ of the speckles of scattered light formed on the CCD chip is given by

$$sd\lambda /a$$
(3)
where $`d`$ is the distance between the diaphragm and the CCD chip, $`\lambda `$ is the wavelength of light, and $`a`$ is the diameter of the diaphragm. By adjusting $`d`$ and/or $`a`$, the size of the speckles relative to the size of a pixel on the CCD chip can be controlled. We use a 2 mm $`\times `$ 2 mm CCD chip with $`128\times 128`$ pixels. Each pixel occupies approximately $`256\mu \text{m}^2`$ with a 100% fill factor.

To maximize the signal-to-noise, we need to collect data from as many speckles of scattered light as possible. However, if the speckles are too small, for example if several speckles occupy a single pixel, optical contrast is lost and the signal for determining a temporal autocorrelation function $`g_2(t)`$ is degraded. Therefore, it is important to adjust the speckle size so that we obtain a satisfactory optical contrast with as much speckles as possible analyzed.

In order to determine the number of speckles per pixel, we adjust the collection optics (i.e. the size of the iris) and measure the spatial intensity autocorrelation function

$$C(p)=\frac{I_iI_{i+p}_i}{I_i^2_i},$$
(4)
where $`I_i`$ is the intensity of scattered light at pixel $`i`$ and $`I_{i+p}`$ is the intensity of light $`p`$ pixels away. The correlation function $`C(p)`$ is determined by averaging over all pixels. Figure 2 shows $`C(p)`$ for three different sets of light collections optics. Note that $`C(p)`$ first decreases and then oscillates about a baseline intensity. We arbitrarily define the position of the first minimum to be the average speckle size. Hence the speckle size is defined as the distance over which the intensity profile remains spatially correlated.

To determine optimal configuration of the collection optics, we use the static scattering pattern from an immobile glass frit. The frit is sufficiently opaque that diffusion approximation for light transport through the medium is valid. The intercept for the intensity correlation function (i.e. $`g_2(0)`$) is expected to lie between 1 and 2. In order to quantitatively evaluate its value, a series of 2000 intercepts was measured. Measurements were performed by self-correlating 2000 images from the glass frit. For each image the mean intensity was calculated as well as the mean squared intensity. The ratio of theses two quantities corresponds to $`g_2(0)`$ for this image. Six series of such measurements were made, each for a different speckle size. Figure 3(a) shows the six distributions of $`g_2(0)`$. The average value $`g_2(0)`$ for each speckle size is an evaluation of the contrast, whereas the spread $`\mathrm{\Delta }g_2(0)`$ indicates how quickly (i.e. how many images must be correlated before) $`g_2(0)`$ converges to its proper statistical value. As such, $`\mathrm{\Delta }g_2(0)`$ quantifies the statistical and instrumental noise. To determine the best compromise between obtaining the highest signal with the lease noise, we plot $`g_2(0)1/\mathrm{\Delta }g(0)`$ as a function of speckle size. Figure 3(b) shows that the best configuration is obtained for a speckle size of 3 pixels in diameter. All subsequent measurements are performed in this configuration. The average intercept $`g_2(0)`$ allows the determination of the normalization factor given for the Siegert relation:

$$g_2(t)1=\beta |g_1(t)|^2$$
(5)
where $`g_1(t)=E(t_0)E(t_0+t)/|E(t_0)|^2`$ is the normalized electric field autocorrelation function. We take $`\beta =0.82\pm 0.05`$ in subsequent measurements and display our results in terms of the squared electric field correlation function $`|g_1(t)|^2`$ using the Siegert relation above.

## IV Calculation algorithm

In this section we explain how the bitmapped image is processed to calculate the autocorrelation function.

### IV.1 Correlation Function Calculation

At time $`t`$, a single CCD frame of scattered light is saved in the computer memory (RAM). We call $`I_i^{raw}(t)`$ the raw intensity of the light incident on pixel $`i`$ for the image taken at time $`t`$, and subtract off an average background $`I_i^{back}`$ for pixel $`i`$ due to stray light and dark counts:

$$I_i(t)=I_i^{raw}(t)I_i^{back}.$$
(6)
We measure $`I_i^{back}`$ by averaging 1000 images for each pixel with the laser turned off. Further refinements taking into account other sources of extraneous signal can be found in cipelletti99 . While such refinements prove useful for DLS in the single scattering limit, we found that such corrections to be insignificant for DWS.

One practical problem encountered in acquiring data is that some images may be skipped by the acquisition card due to synchronization problems. Thus, as a precaution, we always record the time $`t`$ at which each image is obtained.

Applying Eq. (6) to our data, we calculate the correlation function:

$$g_2^{(m)}(t_0,t)=\frac{I_i(t_0)I_i(t+t_0)_i}{I_i(t_0)_iI_i(t+t_0)}_i$$
(7)
Note that the average is performed over *pixels* rather than using the customary method of averaging over the initial time $`t_0`$. That is, the correlation function $`g_2^{(m)}(t_0,t)`$ is calculated by multiplying two frames (one at time $`t_0`$, the other at $`t_0+t`$) together pixel by pixel, summing the results from all pixels, and then dividing by the number of pixels. In general, *one frame suffices to obtain a good average* for $`g_2^{(m)}(t_0,t)`$ because it consists of a large number of pixels, $`16384=128^2`$ in this case, which constitutes several thousand independent samples from which to determine $`g_2^{(m)}(t_0,t)`$.

The correlation function $`g_2^{(m)}(t_0,t)`$ is obtained in this manner is obtained from an *ensemble average* of the sample, with no need for the time averaging method typically used with conventional light detection. When data is taken from an ergodic system, both methods give the same result, as illustrated in Fig. 4 where we compare correlation functions obtained from a CCD camera to one taken with a conventional single detector-single correlator for a sample of polystyrene beads in glycerol obtained undergoing Brownian motion. Significant advantages accrue for the CCD technique, however, when measurements are performed on a non-ergodic sample. Moreover, the CCD method can provide important advantages when measuring correlation functions for non-ergodic systems.

### IV.2 Multiple correlation function calculation

In the preceding section, we discussed how a two-time correlation function $`g_2(t_0,t)`$ can be obtained, $`t_0`$ being the time at which the reference image is acquired and $`t_0+t`$ being a later time at which another image is collected and correlated with the reference. Such a correlation function provides a snapshot of the dynamical state of the system at a particular starting time, namely $`t_0`$, and is possible only because sufficient ensemble averaging is achieved by correlating thousands of pixels between two frames from a CCD camera. The ability to obtain such snapshots is particularly useful for studying systems whose dynamics are time-dependent, i.e. non-stationary systems. Examples include systems which are gelling, coarsening, aging, or undergoing some other non-equilibrium non-stationary process. By measuring $`g_2(t_0,t)`$ as a function of the starting time $`t_0`$, the evolving dynamics of such systems can be followed and studied.

Here, we explain a scheme for following the evolving dynamics of a system by measuring many different correlation functions in parallel, each with a different starting time $`t_0`$. A unique and important feature of this scheme is that new correlation functions can be started at intervals as short as 1/500 s, while still collecting frames for correlation functions started at earlier times. In fact, up to 500 correlation functions can be collected simultaneously, each with a different starting time $`t_0`$, and without any requirement that a previously started correlation function be finished before starting a new one.

The key to the success of the scheme is the realization that in most cases, it is not necessary to measure the correlation function with the same temporal resolution for all delay times. Typically, small time steps are required for short delay times while larger time steps are required for longer delay times. The scheme we employ takes advantage of this fact. It is essentially an interlacing scheme in which successive frames are sometimes used in different correlation functions, as will be explained below. In effect, the scheme represents a compromise between making maximal use of the information available in the measurements for calculating the correlation functions and obtaining and processing the data from the CCD camera frames at a rapid rate.

This compromise is necessary because calculating the correlation function between two frames is time consuming; only one pair of frames can be processed (multiplied and averaged) between the acquisition of two successive images. Image stacking is not practical because it quickly leads to memory overflow. Therefore, we employ an algorithm for calculating many correlation functions, labeled by $`j=1,2,3,\mathrm{},j_{max}`$, simultaneously in real time, that makes only the following manageable demands on the computer:

* only one frame multiplication is performed between the acquisition of successive frames,
* storage of a maximum of $`j_{max}500`$ reference image frames plus the most recently acquired frame in memory, and
* storage of a maximum of $`j_{max}`$ correlation functions in memory, each with on the order of 100 channels.

For a frame rate of 500 frames/s from a $`128\times 128`$ pixel CCD camera, these demands can be met by a computer with a Pentium processor running at 500 MHz and possessing 512 MB of RAM. An important feature of the algorithm is that it requires that only the reference frame for each correlation function, along with most recently acquired frame, be stored in the computer at any moment in time; this keeps the number of frames that must be stored in memory to a manageable number. This also allows all correlation functions to be calculated in real time as the experiment proceeds.

The scheme we use is illustrated graphically in Fig. 5. The vertical axis represents the different correlation functions labeled by $`j=1,2,3,\mathrm{},j_{max}`$. The horizontal axis represents the successive frames that are acquired, labeled $`n=1,2,3,\mathrm{}`$ with each frame being nominally separated in time by the reciprocal $`1/f`$ of the frame rate. The reference image for the $`j^{th}`$ correlation function is given by image number

$$n=n_1^{(j)}=1+(N+1)(j1).$$
(8)
where the subscript “1” indicates that this is the first frame of the $`j^{th}`$ correlation function. Each of these reference images is stored in memory; $`N`$ is typically of order 10. The reference image and the next $`N1`$ images in sequence are correlated with the reference image. The remaining frames used to calculate the $`j^{th}`$ consists of a sequence of selected images as illustrated in Fig. 5. The sequence numbers of the frames used to calculate the $`j^{th}`$ correlation function is summarized by the formula:

$$n=n_i^{(j)}=\{\begin{array}{cc}n_1^{(j)}+i1=(N+1)(j1)+i\hfill & \text{ if }1iN\hfill \\ (j1)+\frac{1}{2}i(i1)\frac{1}{2}N(N1)\hfill & \text{ if }i>N\hfill \end{array}.$$
(9)
Note that the spacing between the points in a correlation function is constant over the first $`N`$ points but increases quadratically at later times. The time interval between two successive correlation functions $`j`$ and $`j+1`$ is approximately $`N/f`$ where $`f`$ is the frame rate. This time interval can be tuned in order to best accommodate the rate of evolution of the system under study. For $`N=1`$ this time interval can be as short as the inverse frame rate.

In Fig. 6 we show a series of correlation functions obtained using the multispeckle CCD camera technique described above. The correlation functions were obtained from a sample consisting of 0.5 wt.% alumina (Al2O3) beads, approximately 160 nm in diameter (AKP30 Sumitomo), that were embedded in a crosslinked elastic network of polybutylacrylate. The first correlation function was started at a time $`t_0=0`$ immediately after a step strain of 5% was applied to the sample. Thus, the correlation functions measure the subsequent random motion of the beads that arises as a response to the stress within the sample. As the stress relaxes over the course of tens of minutes, the random motion of the beads becomes slower. This slowing down is reflected in the increased decorrelation time of correlation functions started at increasingly later times after the initial step strain. In Fig. 7 we plot the decorrelation time $`\tau _{1/2}`$, defined as the time it takes a correlation function to decay to half its initial value, vs. time, referenced to the starting time $`t_0`$. It exhibits a power law behavior over two decades. A detailed discussion of this experiment will be published elsewhere.

## V Fast and slow time scales

The shortest decorrelation time $`\tau _{1/2}`$ measured in the experiments described in the previous section is fairly slow, approximately 0.1 s. In this particular instance, the frame rate of the camera is sufficient to capture all the dynamics of the process in the correlation function. This is not always the case, however. Many materials, including polymer networks, colloidal gels, and glassy materials possess two characteristic ranges of relaxation times, where only the slower of the two is sufficiently slow for the multispeckle CCD camera method to work. In these cases, the faster relaxation processes are not measured.

To illustrate such a process, we use the multispeckle technique to follow the gelation that occurs following a quench from $`60^{}`$C to $`20^{}`$C in a 1% gelatin solution containing 1% latex spheres with a radius of 85 nm. The correlation functions are shown in Fig. 8. As the gelatin gels, the motion of the strands of gelatin and the motion of the spheres embedded within it becomes increasingly restricted. As the gelation proceeds, the height of the plateau and the $`y`$-intercept of the correlation functions increases. This occurs because there is a fast relaxation process that occurs outside of the time window of the camera. That is, there is an initial rapid drop of $`g_2(t)`$ that occurs on a time scale that is faster than the inverse frame rate of the camera ($`2`$ ms). Below, we explain how this obstacle can be overcome such that the lower bound of the time window is reduced from $`10^2`$ to $`10^8`$ s.

What we would like to be able to do is to *simultaneously* capture the fast and slow dynamics for a time-dependent non-stationary process. Below, we show that this can be accomplished, provided at strict set of criteria are met. Fortunately, it turns out that these criteria are met by most of the systems of interest so that the criteria are not very restrictive in practical situations. As we will show below, we can use a classical photomultiplier-based (fast) correlator, which spans a delay-time range of $`10^8`$ s to about $`1`$ s, to capture the fast dynamics and the multispeckle (slow) correlator to capture the slow dynamics. The key is to properly ensemble average the scattered light signal entering the photomultiplier used with the fast correlator without disturbing the multispeckle signal that enters the slow correlator.

To provide an ensemble-averaged signal to the fast correlator, we exploit a technique recently developed by Scheffold et al. romer00 , which enables simple and efficient spatial averaging without requiring any motion of the sample or the detector. This technique is known as the two-cell technique and has been described both theoretically skipetrov96  and experimentally romer00 .

### V.1 Two-cell technique

The idea of the two-cell technique is to optically couple the studied sample to an ergodic reference system. The light is successively scattered by the sample cell and then by the reference cell before entering the photomultiplier. This can be implemented simply by illuminating two cells in series, and detecting the light that passes through and is scattered by both samples in the two cells. In that case, the non-ergodic signal emerging from the sample cell is made ergodic by its scattering through the reference cell. The second cell can thus be regarded as performing a slow spatial averaging of the sample signal. Under carefully specified optical conditions, the correlation function, $`g_1^M(t),`$ measured after passing through the two cells is simply the product of the correlation function of each cell:

$$g_1^M(t)=g_1^S(t)\times g_1^R(t),$$
(10)
where the superscripts $`R`$ and $`S`$ stand for the sample and the reference correlation functions, respectively. The conditions under which the two-cell technique can be employed have been carefully considered in Ref. skipetrov96 . They find that the primary condition required is that each photon traverses the space between the sample and reference cells only once. This means that for a configuration where a significative amount of light is lost between the two cells, as was done in Ref. romer00 , this condition is satisfied when $`l_R^{}/L_Rl_S^{}/L_S`$, where $`l_i^{}`$ and $`L_i`$ are, respectively, the photon transport mean free path and the cell thickness for the reference ($`i=R`$) and the sample ($`i=S`$) cells.

### V.2 Experimental set up

Figure 9 illustrates how we incorporate the two-cell technique into the multispeckle correlator so that we can simultaneously capture the fast and slow dynamics from a slowly evolving non-ergodic system. We arrange the optics for the multispeckle (slow) correlator just as before so that surface from which the multiply scattered light exits the sample cell is imaged with a 1:1 ratio onto a pinhole behind which the CCD camera is placed. For the fast correlator, a non-polarizing beam splitter is placed a short distance in front of the pinhole so that half of the light is diverted onto the reference cell. The reference cell positioned so that the exit plane of the sample cell is imaged onto the reference cell’s front surface. An single mode optical fiber with a collimator collects the scattered light exiting the reference cell at an angle of approximately $`15^{}`$ and sends it to a photomultiplier whose output is connected to a multiple-sample-time correlator.

The reference cell has a thickness of 5 mm and contains a glycerol solution of 0.5 $`\mu `$m-radius latex spheres with a volume fraction of $`1\%`$. The concentration of spheres is chosen so that there is just enough multiple scattering in the reference cell to ensure that essentially no light passes through the sample without being scattered at least several times but with $`l_R^{}/L_R`$ as big as possible. However, this must be done so that the decorrelation time of the reference correlation function is somewhat larger, typically by a factor of about 3 to 10, than the inverse frame rate of the CCD camera. By adjusting the concentration of spheres and the viscosity of the solution with glycerol, these conditions are readily achieved.

Before a measurement can be made using the two-cell technique, it is necessary to measure the correlation function of the reference cell in the absence of the sample cell and to ensure that its correlation time is greater than the inverse frame rate. Data from such a measurement is shown in Fig. 10, from which one can see that the decay time of $`g_1^R(t)`$ is about 0.1 s. The correlation function $`g_1^M(t)`$ is then measured with the sample cell in place, from which the sample correlation function $`g_1^S(t)`$ is obtained using Eq. (10). Figure 10 shows results from a typical set of measurements, including the reference cell correlation function $`g_1^R(t)`$, the measured two-cell correlation function $`g_1^M(t)`$, and the true correlation function of the sample $`g_1^S(t)`$ obtained from Eq. (10). Note that $`g_1^S(t)`$ is obtained accurately from the first time measured by the correlator, about $`10^7`$ s, out to about 30 ms, or approximately the decorrelation time of $`g_1^R(t)`$. For longer delay times, the reference correlation function $`g_1^R(t)`$ falls to zero making recovery of $`g_1^S(t)`$ using Eq. (10) impossible. Note that the sample correlation function $`g_1^S(t)`$ determined in this way has a significant time interval of overlap, from about 2 ms to about 30 ms, with the correlation function measured by the multispeckle technique. Thus, it should be possible to superimpose the correlation functions from the two techniques to obtain a single correlation function covering the entire range of time scales from 10 ns to hundreds of seconds or longer. It is important to note, however, that in order from the signal to be ergodic the correlation function obtained from the two-cell technique to be properly ensemble averaged, data must be collected over a time that is about $`10^3`$ times larger than the decay time of the reference correlation function. In our case the acquisition time is 100 s.

### V.3 Overlapping data from the two techniques

In order to overlay the correlation functions measured by the two techniques, they must be properly normalized using the respective Siegert coefficients $`\beta `$ for each measurement (see Eq. 5). For the multispeckle CCD camera method, $`\beta =0.8`$, as discussed in Section III.1. For the double cell technique, $`\beta `$ is determined by the value of the intercept from the correlation function of the reference cell alone. For the reference cell, we obtain $`\beta =0.95`$.

Figure 11 shows correlation functions from the two techniques for the gelation experiment described in Section V. The acquisition time for each double cell measurement was 100 s. Excellent agreement is obtained over the range where the two techniques overlap (2 ms-30 ms). Note that the overlap occurs in a region where the slope of the correlation function is sizeable, making the agreement even more impressive.

Finally, it us useful to compare the results for different methods of measuring a correlation function. In Fig. 12, we plot three correlation functions: (1) one obtained using the multispeckle technique with a window of decay times of 2 ms–10 s. Averaging was performed in real time (up to the longest decay time or $`20`$ s); (2) another obtained using the double cell technique over the decay-time window of $`10^7`$s–30 ms and averaged for 100 s; (3) and a third obtained using a classical set-up with a photomultiplier tube and single correlator over the decay-time window of $`10^75`$ s and averaged for 3000 s. The measurements using the double cell and the multispeckle methods were done simultaneously; the measurement with the photomultiplier and single correlator was performed just after the other two, simply by removing the reference cell. The agreement is very good over the time window of $`10^7`$-$`10^1`$ s. There are deviations, however, between the single correlator and the multispeckle techniques at longer delay times. These deviations occur because data obtained using the single correlator technique must be collected over a fairly long data-acquisition time, 100 s in this case, in order to obtain satisfactory signal-to-noise. In cases where the dynamics of the sample evolves with time, the multispeckle technique has clear advantages.

The multi-speckle technique can measure correlation functions in real time for delay-times larger than 2 ms. When coupled simultaneously (in the same set-up) with the two-cell technique, the combined method can measure correlation functions over an extremely broad range of delay times. As a demonstration, we show in Fig. 13 a correlation function measured for a gel of silica beads (0.2 $`\mu m`$ in diameter) at a volume fraction of 52% with an Ionic force adjusted to $`1.4`$ $`10^{}2`$M with KNO3. correlation function spans a range of delay times from 10$`{}_{}{}^{8}s`$ to greater than 103 s, or more than 11 decades in time. Such a measurement is not only remarkable, it is extremely useful for probing glassy systems which can exhibit dynamics over a vary broad range of delay times.

Lastly, there are a couple of practical constraints that bear mentioning. First, the two techniques can be usefully coupled only if the typical time scale over which the system evolves is much longer than the averaging time required for the two cell technique. As a practical matter, the required averaging time is set by the longest delay time one wishes to measure. That time is determined primarily by the desire to have some range of delay times where the two-cell and multi-speckle techniques overlap. The limitation here is the speed of the camera that can be used for the multi-speckle technique. Second, the simultaneous use of the two techniques leads to a certain amount of light loss. Thus, a powerful laser ($`1`$ W) is recommended.


format: text

text:
# Untitled Document

SUPERFLUID PHASES OF TRIPLET PAIRING

J. W. Clark

McDonnell Center for the Space Sciences

and Department of Physics, Washington University

St. Louis, Missouri 63130, USA

V. A. Khodel and M. V. Zverev

Russian Research Center Kurchatov Institute

Moscow 123182, Russia

1. INTRODUCTION

Triplet pairing is exemplified physically in at least two strongly interacting quantum many-fermion systems:

(i) In laboratories on earth: Superfluid phases are observed when liquid 3He is cooled to mK temperatures (triplet pairing of spin–1/2 atoms) .

(ii) In the quantum fluid interior of a neutron star: It is generally believed that the neutrons pair-condense at temperatures below $`10^9`$ K (triplet pairing of spin-1/2 nucleons) .

There are important differences in these realizations. In dense neutron matter, $`{}_{}{}^{1}S_{0}^{}`$ pairing, which dominates in the inner crustal region of the star, is quenched by the strong short-range repulsion in the nucleon-nucleon ($`NN`$) interaction. At $`NN`$ collision energies corresponding to the relevant density regime, the most attractive phase shift is found in the $`{}_{}{}^{3}P_{2}^{}`$ channel, indicating that this state will dominate the pairing. The $`{}_{}{}^{3}P_{0}^{}`$ phase shift is only weakly attractive, while the $`{}_{}{}^{3}P_{1}^{}`$ phase shift is repulsive; these channels are therefore expected to have only minor effects on the problem. The state dependence of the spin-orbit component of the $`NN`$ interaction is responsible for this pre-eminence of the $`{}_{}{}^{3}P_{2}^{}`$ state. A complication is introduced by the tensor component of the interaction, which implies a coupling of the $`{}_{}{}^{3}P_{2}^{}`$ and $`{}_{}{}^{3}F_{2}^{}`$ channels.

On the other hand, in liquid $`{}_{}{}^{3}He`$, the atom-atom interaction is spherically symmetric and spin-independent to a very good approximation. Tensor and spin-orbit components are effectively absent. Thus there is no coupling to higher $`L`$ waves, but the three $`J`$ states for $`L=1`$, $`S=1`$, namely $`{}_{}{}^{3}P_{2}^{}`$, $`{}_{}{}^{3}P_{1}^{}`$, and $`{}_{}{}^{3}P_{0}^{}`$ now enter on an equal footing.

On balance, the neutron matter – in spite of the complication of tensor coupling – is simpler than liquid $`{}_{}{}^{3}He`$, since fewer magnetic substates enter the problem. Therefore we focus first on neutron matter as a “training ground” for the more challenging case of liquid 3He.

Quantitative, microscopic prediction of the superfluid phases and phase diagram of systems manifesting triplet pairing (or pairing in orbital angular momentum states beyond $`S`$-wave), has hitherto eluded the best efforts of theorists. For one thing, the popular Ginzburg-Landau approach is restricted to temperatures near the critical temperature $`T_c`$. Secondly, iterative procedures, routinely employed to calculate the energy gap in weakly interacting systems with $`S`$-wave pairing, are afflicted with slow convergence and uncertain accuracy when applied to the many coupled, singular nonlinear integral equations that come into play for pairing in higher angular momentum states. The limitations of standard iterative procedures become most apparent when one attempts to construct the superfluid phase diagram of the system, which is controlled by tiny energy splittings between different solutions of the BCS pairing problem.

These difficulties can now be overcome by means of a recently introduced separation method \[3–8\] for solving the BCS gap equations associated with pairing states of arbitrary angular-momentum content. In this paper, we implement this method to identify and classify the spectrum of phases of the $`{}_{}{}^{3}P_{2}^{}`$$`{}_{}{}^{3}F_{2}^{}`$ pairing model. Nominally, the model refers to superfluid neutron matter at densities of order twice that found in heavy nuclei. However, the results obtained are universal in the sense that the angular properties of the solutions corresponding to the existing phases are independent of the detailed nature of the system under study, and especially its pairing interaction. Indeed, these results can be derived analytically and depend only on two energy scales, one governing the overall strength of the pairing effect and the other measuring the strength of the coupling between $`P`$ and $`F`$ states. There exist the usual solutions involving a single value of the magnetic quantum number $`M`$ (and its negative) – obviously three in number corresponding to $`M=0,\pm 1,\pm 2`$. In addition we find that there exist ten real multicomponent solutions. Five of these have angle-dependent order parameters with nodes (and therefore are of relatively high energy), while the other five are nodeless. In contrast to the case of superfluid 3He, transitions occur between phases with nodeless order parameters. Results for the temperature dependence of the competition between the various phases have been obtained. In principle, and presumably also in practice, the same approach can be used to complete the catalog of superfluid phases of liquid 3He.

2. BCS FORMALISM FOR ANY-CHANNEL PAIRING

Adopting a spinor representation, the $`2\times 2`$ gap matrix for the general BCS problem of pairing in a state of triplet spin ($`S=1`$) and triplet isospin ($`T=1`$) has the expansion

$$\mathrm{\Delta }_{\alpha \beta }(𝐩)=\underset{J,L,M}{}\mathrm{\Delta }_L^{JM}(p)\left(G_{LJ}^M(𝐧)\right)_{\alpha \beta }$$
$`(1)`$
in terms of the spin-angle matrices

$$\left(G_{LJ}^M(𝐧)\right)_{\alpha \beta }=\underset{M_SM_L}{}C_{\frac{1}{2}\frac{1}{2}\alpha \beta }^{1M_S}C_{1LM_SM_L}^{JM}Y_{LM_L}(𝐧).$$
$`(2)`$
The coupled partial-wave gap components $`\mathrm{\Delta }_L^{JM}(p)`$ solve the set of gap equations

$$\begin{array}{ccc}\hfill \mathrm{\Delta }_L^{JM}(p)& =\underset{L^{}L_1J_1M_1}{}(1)^{1+\frac{LL^{}}{2}}p|V_{LL^{}}^J|p_1S_{L^{}L_1}^{JMJ_1M_1}(𝐧_1)\hfill & \\ & \times \frac{\mathrm{tanh}\left(E(𝐩_1)/2T\right)}{2E(𝐩_1)}\mathrm{\Delta }_{L_1}^{J_1M_1}(p_1)p_1^2dp_1d𝐧_1.\hfill & (3)\hfill \end{array}$$

The quasiparticle energy $`E(𝐩)=\left[\xi ^2(p)+D^2(𝐩)\right]^{1/2}`$ is constructed from the gap components $`\mathrm{\Delta }_L^{JM}(p)`$ through

$$D^2(𝐩)=\underset{LJML_1J_1M_1}{}\left(\mathrm{\Delta }_L^{JM}(p)\right)^{}\mathrm{\Delta }_{L_1}^{J_1M_1}(p)S_{LL_1}^{JMJ_1M_1}(𝐧)$$
$`(4)`$
together with the single-particle spectrum $`\xi (p)`$ of the normal Fermi liquid, often parametrized in terms of an effective mass $`M^{}`$. Angular dependence is introduced into the quasiparticle energy $`E(𝐩)`$ via the spin trace

$$S_{LL_1}^{JMJ_1M_1}(𝐧)=\mathrm{Tr}\left[\left(G_L^{JM}(𝐧)\right)^{}G_{L_1}^{J_1M_1}(𝐧)\right],$$
$`(5)`$
which obviously complicates explicit solution of the system of gap equations.

The pairing matrix elements $`p|V_{LL^{}}^J|p_1`$ are generated by the spin-angle expansion

$$V(𝐩,𝐩_1)=\underset{LL^{}JM}{}(1)^{\frac{LL^{}}{2}}p|V_{LL^{}}^J|p_1G_{LJ}^M(𝐧)\left(G_{L^{}J}^M(𝐧_1)\right)^{}$$
$`(6)`$
of the totality of vertex diagrams irreducible in the particle-particle channel.

A close inspection reveals that complete solution of this set of equations presents awesome difficulties: note, in particular, the coupling to the generic angular momentum labels $`L_1`$, $`J_1`$, and $`M_1`$ quantum numbers via the squared-gap quantity $`D^2(𝐩)`$ appearing in the quasiparticle energy $`E(𝐩)`$. Practical approximate solution (which might still be extremely accurate) will require a series of justifiable simplifications.

To wit, in dealing with neutron matter, the free-space two-body interaction has the salient feature that the components of the central forces nearly cancel each other, as reflected in the behavior of the experimental $`P`$-scattering phases. It is assumed that this feature is preserved by the effective interaction inside neutron matter. We shall later comment on the veracity of this assumption. If it holds, the pivotal role of the spin-orbit force in promoting the $`{}_{}{}^{3}P_{2}^{}`$ pairing channel then implies that contributions to triplet pairing from “nondiagonal” terms with $`L^{},L_11`$ or $`J_12`$ on the right-hand side of the set of gap equations can be evaluated within perturbation theory, in terms of the set of principal gap amplitudes $`\mathrm{\Delta }_1^{2M}(p)`$, with $`M=0,\pm 1,\pm 2`$.

In fact, if we choose to invoke time-reversal invariance, the problem may be treated in terms of only three complex functions, namely $`\mathrm{\Delta }_1^{2M}(p)`$ with $`M=0,1,2`$.

Another simplifying feature is the existence of a small parameter across the range of interesting pairing problems. When dealing with anisotropic gaps as may arise in pairing states beyond the $`S`$-wave, it is helpful to define the energy gap parameter $`\mathrm{\Delta }_F`$ as the square root of the angle average of $`D^2(𝐩)`$, evaluated at the Fermi surface. The ratio $`d_F=\mathrm{\Delta }_F/ϵ_F`$ of the gap to the Fermi energy provides a small parameter for BCS pairing theory. For $`{}_{}{}^{3}P_{2}^{}`$ (or $`{}_{}{}^{3}P_{2}^{}`$$`{}_{}{}^{3}F_{2}^{}`$) pairing in neutron matter, a leading order approximation in $`d_F`$ is good to about one part in $`10^5`$ or $`10^6`$. A similar or higher accuracy is to be expected in liquid 3He.

3. THE SEPARATION METHOD

The separation approach developed in Refs. \[3-8\] facilitates essential further simplifications (linked in part to the existence of a (very) small parameter). In this approach, any given pairing matrix element is expressed identically as a separable part, plus a remainder that vanishes when either momentum argument is on the Fermi surface. We write

$$p|V_{LL^{}}^J|p_1=v_{LL^{}}^J\varphi _{LL^{}}^J(p)\varphi _{LL^{}}^J(p_1)+W_{LL^{}}^J(p,p_1)$$
$`(7)`$
and verify that the choices $`v_{LL^{}}^Jp_F|V_{LL^{}}^J|p_F`$ and $`\varphi _{LL^{}}^J(p)p|V_{LL^{}}^J|p_F/v_{LL^{}}^J`$ produce the desired behavior when $`p`$ or $`p_1`$ hits $`p_F`$. Substitution of this identity into the set of gap equations, followed by simple manipulation and argumentation, establishes a decomposition of the form

$$\mathrm{\Delta }_L^{JM}(p)=D_L^{JM}\chi _{LL^{}}^J(p)(|M|=0,1,2)$$
$`(8)`$
for the general gap component and provides separate equations for the shape factor $`\chi _{LL^{}}^J(p)`$ (normalized to unity at $`p=p_F`$) and the numerical amplitudes $`D_L^{JM}`$ that determine the angle-dependence of the gap.

The shape factor $`\chi _{LL^{}}^J(p)`$ is given by a nonsingular integral equation involving a set of kernels proportional to the residual interaction $`W_{LL^{}}^J(p,p_1)`$ (see Ref.  for details). The vanishing of $`W`$ at the Fermi surface removes the logarithmic singularity characteristic of BCS theory, which is banished to the set of equations for the amplitudes $`D_L^{JM}`$, which are coupled nonlinear equations for a set of numbers.

To very high accuracy, the integral equation for $`\chi _{LL^{}}^J`$ is linear and independent of the amplitudes $`D_L^{JM}`$. For all practical purposes (e.g., appealing to extremely rapid convergence of small-parameter expansions), we are free to make the replacements

$$E(𝐩)|\xi (𝐩)|,\mathrm{tanh}(E(𝐩)/2T)1$$
$`(9)`$
in any integral involving the residual interaction $`W`$ as a factor. The shape factors $`\chi _{LL^{}}^J(p)`$ are therefore determined by the interaction $`W`$ and may be calculated independently of the $`D_L^{JM}`$ by a matrix inversion routine. The nonlinear and singular aspects of the problem reside entirely in the set of equations for the amplitudes $`D_L^{JM}`$, which may also be solved numerically by standard methods (e.g. Newton-type algorithms). In fact, we have even been able to solve scaled versions of the $`D_L^{JM}`$ equations analytically in interesting cases, as will emerge below.

In view of these crucial simplifications revealed and expedited by the separation transformation, it is seen that

(i) For given $`L`$, $`L^{}`$, and $`J`$, the $`\chi `$ function is universal, i.e., independent of magnetic quantum number $`M`$ and temperature $`T`$.

(ii) The factorization

$$\mathrm{\Delta }_1^{2M}(p)=D_1^{2M}\chi (p)(|M|=0,1,2)$$
$`(10)`$
may be asserted for the principal gap components, where $`\chi (p)\chi _{11}^{21}(p)`$.

(iii) Explication of the phase diagram of dense superfluid neutron matter reduces to determination of the three amplitudes $`D_1^{2M}`$, since the character of the phase diagram itself is independent of the shape factor $`\chi (p)`$.

4. PERTURBATIVE TREATMENT

Now let us focus on the nondiagonal contributions to the right-hand side of the gap equations (3), where nondiagonal means one or more of the conditions $`L^{}=1`$, $`L_1=1`$, $`J_1=2`$ is not met. Two nondiagonal contributions are found to be of leading importance: The first contains the integral of the product $`V_{31}^2S_{31}^{2M2M_1}\mathrm{\Delta }_1^{2M_1}`$ while the second contains the integral of the product $`V_{11}^2S_{13}^{2M2M_1}\mathrm{\Delta }_3^{2M_1}`$. There is a single small factor in each product, namely $`V_{31}`$ in the first instance and $`\mathrm{\Delta }_3^{2m_1}`$. All other nondiagonal contributions involve two or more small factors. Retaining the leading nondiagonal pieces and ignoring all the others, we have what is called the $`{}_{}{}^{3}P_{2}^{}`$$`{}_{}{}^{3}F_{2}^{}`$ pairing problem, which is generally regarded as a satisfactory model of superfluid neutron matter.

Rapid convergence of the pertinent nondiagonal integrals is instrumental to success of the perturbation approach used to treat the nondiagonal effects. To affirm this behavior, we observe that the dominant contributions to these integrals come from the vicinity of the Fermi surface. Then, if $`E(𝐩)`$ is significantly larger than the energy gap value (as it will be for large $`p`$), $`E(𝐩)`$ and $`|\xi (𝐩)|`$ will practically coincide, and consequently the angular integration yields zero.

We conclude that in evaluating the nondiagonal contributions it suffices to know the “minor” gap components $`\mathrm{\Delta }_3^{2M}(p)`$ at $`p=p_F`$. These quantities will in fact be determined in terms of the coefficients $`D_1^{2M}`$. To excellent approximation, we may retain only the dominant contribution to the right-hand side of the gap equations (3) that diverges like $`\mathrm{ln}d_F`$ as $`\mathrm{\Delta }_F0`$. We are then led to the desired connection

$$\mathrm{\Delta }_3^{2M}(p=p_F)=\eta D_1^{2M},$$
$`(11)`$
where $`\eta =p_F|V_{13}^2|p_F/v_F`$ and $`v_Fp_F|V_{11}^2|p_F`$. Similar relations may be obtained for other minor gap components, notably $`\mathrm{\Delta }_1^{00}`$ and $`\mathrm{\Delta }_1^{1M}`$.

What can we say about the parameter $`\eta `$, which represents the coupling to the “nondiagonal” states? If perturbation theory is to be valid, it should be comfortably small. For pairing matrix corresponding directly to in-vacuum neutron-neutron interaction, $`\eta `$ depends smoothly on density, varying around 0.3 in the interval $`\rho _0<\rho <2\rho _0`$. Medium modification of the spin-orbit force is probably not important. Due to the relativistic origin of this component, it should not be much affected by polarization or correlation corrections. This claim is consistent with empirical analyses of the spin-orbit splitting in finite nuclei.

On the other hand, medium modification of the tensor force may be more significant, especially as one approaches the density at which pion condensation occurs. Thus, the parameter $`\eta `$ is somewhat uncertain, although it is still expected to be rather small.

To sketch out the superfluid phase diagram of the system, we need to determine the key amplitudes $`D_1^{2M}`$ ($`M=0,1,2`$) to leading perturbative order in the coupling $`\eta `$. Exploiting the linear connection (11) between the $`D_1^{2M}`$ amplitudes and the minor-component values $`\mathrm{\Delta }_3^{2M}(p_F)`$, simple manipulations applied to the coupled gap equations (1) at $`p=p_F`$ yield

$$D_1^{2M}+v_F\underset{M_1}{}D_1^{2M}\varphi (p)\frac{\mathrm{tanh}\left(E_0(𝐩)/2T\right)}{2E_0(𝐩)}S_{11}^{2M2M_1}(𝐧)\chi (p)p^2𝑑p𝑑𝐧=\eta v_Fr_M$$
$`(12)`$
with $`\varphi (p)p|V_{11}^2|p_F/v_F`$, $`E_0(𝐩)E(𝐩;\eta =0)`$, and

$$r_M=\underset{M_1}{}D_1^{2M_1}\left[S_{31}^{2M2M_1}(𝐧)+S_{13}^{2M2M_1}(𝐧)\right]\frac{\mathrm{tanh}\left(E_0(𝐩)/2T\right)}{2E_0(𝐩)}p^2𝑑p𝑑𝐧.$$
$`(13)`$
We restrict the search for solutions of these equations to those with real coefficients $`D_1^{2M}`$, reasoning that solutions with complex $`D`$ amplitudes will lie at energies high enough to make them physically irrelevant.

Inserting the explicit form of $`S_{11}^{2M2M_1}(𝐧)`$ from Eq. (5) into Eq. (12), we derive a system of three equations for the gap value $`\mathrm{\Delta }_F`$ and the two ratios

$$\lambda _1=D_1^{21}/D_1^{20}\sqrt{6}\mathrm{and}\lambda _2=D_1^{22}/D_1^{20}\sqrt{6},$$
$`(14)`$
which serve generally to determine the angular dependence of the gap (or alternatively its composition with respect to the magnetic quantum number $`M`$). In a notation and form compatible with the earlier treatment of the pure $`{}_{}{}^{3}P_{2}^{}`$ problem , these basic equations read

$$\begin{array}{ccc}\hfill \lambda _2+v_F\left[\lambda _2(J_0+J_5)\lambda _1J_1J_3\right]& =\eta v_Fr_2,\hfill & (15a)\hfill \\ \hfill \lambda _1+v_F\left[(\lambda _2+1)J_1+\lambda _1(J_0+4J_5+2J_3)/4\right]& =\eta v_Fr_1,\hfill & (15b)\hfill \\ \hfill 1+v_F\left[(\lambda _2J_3+\lambda _1J_1)/3+J_5\right]& =\eta v_Fr_0,\hfill & (15c)\hfill \end{array}$$

with

$$J_i=f_i(\theta ,\phi )\varphi (p)\frac{\mathrm{tanh}\left(E_0(𝐩)/2T\right)}{2E_0(𝐩)}\chi (p)\frac{p^2dpd𝐧}{4\pi },$$
$`(16)`$
$$f_0=13z^2,f_1=3xz/2,f_3=3(2x^2+z^21)/2,f_5=(1+3z^2)/2,$$
$`(17)`$
$$z=\mathrm{cos}\theta ,x=\mathrm{sin}\theta \mathrm{cos}\phi ,y=\mathrm{sin}\theta \mathrm{sin}\phi .$$
$`(18)`$
Upon setting $`\eta =0`$, Eqs. (15a)–(15c) collapse to Eqs. (12)–(14) of Ref. , which refer to the case $`\kappa _1=\kappa _2=0`$ and were solved analytically in Refs. .

If we had not used the separation method, we would have been faced at this point with a system of coupled singular nonlinear integral equations for a set of functions, not numbers. The great advantage gained is transparent.

5. SEARCHING FOR MULTICOMPONENT SOLUTIONS

The basic equations (15a)–(15c) possess the familiar one-component solutions corresponding to definite $`|M|`$, i.e., $`|M|=0,1,2`$. In addition, there exist multicomponent solutions whose structure and spectrum we seek to establish, via a two-step process. We first note that $`J_5`$, which is the only $`J_i`$ integral containing a principal term going like $`\mathrm{ln}(ϵ_F/\mathrm{\Delta }_F)`$, is responsible for the gap magnitude $`\mathrm{\Delta }_F`$. On the other hand, $`J_5`$ is irrelevant to the phase structure.

Step 1. Thus, we begin by eliminating $`J_5`$ from the first pair (15a)–(15b) of basic equations and reduce the number of $`J_i`$ integrals in each of the pair to two:

$$\begin{array}{ccc}\hfill (\lambda _2+1)[3\lambda _1(\lambda _2+1)J_02(\lambda _1^22\lambda _2^2+6)J_1]& =\eta B_1,\hfill & (19a)\hfill \\ \hfill (\lambda _2+1)[(\lambda _1^24\lambda _2)J_1+\lambda _1(\lambda _2+1)J_3]& =\eta B_2,\hfill & (19b)\hfill \end{array}$$

$$\begin{array}{ccc}\hfill B_1& =2\lambda _1(2\lambda _2+3)r_24(\lambda _2^23)r_16\lambda _1(\lambda _2+2)r_0,\hfill & \\ \hfill B_2& =\lambda _1r_2+4\lambda _2r_13\lambda _1\lambda _2r_0.\hfill & (20)\hfill \end{array}$$

Step 2. Assuming $`\lambda _21`$, we perform the rotation

$$(x,z)=(t\mathrm{sin}\vartheta +u\mathrm{cos}\vartheta ,t\mathrm{cos}\vartheta +u\mathrm{sin}\vartheta )$$
$`(21)`$
and choose the angle $`\vartheta `$ so as to remove the integral $`J_1`$ from the “sanitized” first pair of basic equations, (19a)–(19b). Specifically, if $`\zeta =\mathrm{tan}\vartheta `$ is chosen to obey the algebraic equation

$$\lambda _1\zeta ^2(\lambda _23)\zeta \lambda _1=0,$$
$`(22)`$
then substitution of the transformed integrals $`J_i`$ converts the key pair of equations (19a)–(19b) into

$$\begin{array}{ccc}\hfill (\lambda _2+1)[A_1J_0+A_2J_3]& =\eta B_1,\hfill & (23a)\hfill \\ \hfill (\lambda _2+1)[A_1J_0+A_2J_3]& =2\eta B_2,\hfill & (23b)\hfill \end{array}$$

with

$$\begin{array}{ccc}\hfill A_1& =\frac{3}{2}\lambda _1(1+\lambda _2)(2\zeta ^2)\frac{3}{2}(\lambda _1^22\lambda _2^2+6)\zeta ,\hfill & \\ \hfill A_2& =3\lambda _1(1+\lambda _2)\zeta ^2(\lambda _1^22\lambda _2^2+6)\zeta .\hfill & (24)\hfill \end{array}$$

We observe that the left-hand members of Eqs. (23a) and (23b) are identical! It is in fact just this coincidence that leads to the striking universalities of the pure $`{}_{}{}^{3}P_{2}^{}`$ pairing problem first revealed in Ref. . Of course, in this problem, which corresponds to $`\eta =0`$, we have the special circumstance that the right-hand side of each equation vanishes identically, so that Eqs. (23a) and (23b) coincide and yield only the single constraint

$$(\lambda _2+1)[A_1J_0+A_2J_3]=0$$
$`(25)`$
on the parameters $`\lambda _1`$ and $`\lambda _2`$.

Figure 1. Solutions for the coefficient ratios $`\lambda _1`$ and $`\lambda _2`$ in the case $`\eta 0`$, corresponding to the uncoupled $`{}_{}{}^{3}P_{2}^{}`$ pairing problem solved in Refs. . The particular (point) solution ($`\lambda _1=0`$, $`\lambda _2=3`$) is indicated by the solid dot and the degenerate solutions by the solid curves.

It ensues that the solutions of the pure $`{}_{}{}^{3}P_{2}^{}`$ pairing problem display remarkable universalities as expressed in two kinds of degeneracies, independently of temperature, density, and details of the in-medium interaction:

$``$Energetic degeneracy. There exists an upper group of states, degenerate in energy, whose angle-dependent order parameters have nodes, and a lower group with nodeless order parameters. Relative to the absolute pairing energy, the splitting between upper and lower states is small (of order 2% in the neutron-matter problem).

Figure 2. Multicomponent solutions of the $`{}_{}{}^{3}P_{2}^{}`$ pairing problem, defined by the coefficient ratios $`\lambda _1`$ and $`\lambda _2`$. Solutions with nodeless (respectively, node-bearing) order parameters in the case of pure $`{}_{}{}^{3}P_{2}^{}`$ pairing are indicated by open (respectively, filled) circles.

$``$Parametric degeneracy. The multicomponent solutions, which satisfy the spectral condition

$$(\lambda _1^2+22\lambda _2)(\lambda _1^22\lambda _2^26\lambda _2)=0,$$
$`(26)`$
display a parametric degeneracy with respect to the parameters $`\lambda _1`$ and $`\lambda _2`$. Thus, as seen in Fig. 1 (which was constructed analytically), the multicomponent solutions generally define curves rather than points in the $`(\lambda _1,\lambda _2)`$ plane.

6. BREAKING THE PARAMETRIC DEGENERACY

The strong parametric degeneracy intrinsic to the uncoupled $`{}_{}{}^{3}P_{2}^{}`$ pairing problem is lifted in the case of $`{}_{}{}^{3}P_{2}^{}`$$`{}_{}{}^{3}F_{2}^{}`$ pairing, i.e., at $`\eta 0`$. With the coupling turned on, the true solutions of the problem are represented by a set of isolated points in the $`(\lambda _1,\lambda _2)`$ plane.

From the analytical standpoint, the underlying “mechanism” runs as follows. Upon equating the right-hand members of the double-sanitized pair of basic equations (23a)–(23b) in the small–$`|\eta |`$ limit ($`|\eta |`$ infinitesimal), one obtains an additional relation between the parameters $`\lambda _1(\eta =0)`$ and $`\lambda _2(\eta =0)`$:

$$\lambda _1r_2(\lambda _23)r_13\lambda _1r_0=0.$$
$`(27)`$
This relation supplements the spectral condition (26) nontrivially and removes the parametric degeneracy. The system formed by (27) and (26) is solved by applying the same rotation in $`xz`$ coordinates as introduced previously. After some algebra, one may then obtain the full set of solutions of the coupled-channel $`{}_{}{}^{3}P_{2}^{}`$$`{}_{}{}^{3}F_{2}^{}`$ pairing problem. A salient feature of these solutions, made explicit via the separation scheme, is their virtually complete independence of the temperature $`T`$.

We summarize the results of this analysis with a catalog of the possible solutions for $`{}_{}{}^{3}P_{2}^{}`$$`{}_{}{}^{3}F_{2}^{}`$ pairing described in the BCS framework. One simplification, already evident in Fig. 1 but also quite general, is that the relevant pairing energies are independent of the sign of $`\lambda _1`$, so that we only need to consider $`\lambda _1>0`$. The modified picture in the $`(\lambda _1,\lambda _2)`$ plane is displayed in Fig. 2.

First of all, there remain the three well-known single-component solutions with $`|M|=0`$, 1, or 2. Beyond these, the collection of unitary solutions of the $`{}_{}{}^{3}P_{2}^{}`$$`{}_{}{}^{3}F_{2}^{}`$ pairing problem contains ten multicomponent solutions, corresponding to more complicated superfluid phases.

Five of these additional solutions, denoted $`O_k`$ ($`k=1,\mathrm{},5`$), have nodeless order parameters and include:

(a) Two two-component solutions $`O_{\pm 3}`$, identical to those found in the pure $`{}_{}{}^{3}P_{2}^{}`$ pairing problem, with $`\lambda _1=0`$ and $`\lambda _2=\pm 3`$.

(b) Three three-component solutions:

$``$Two of them, $`O_1`$ and $`O_4`$, are associated with the upper branch of $`\lambda _1^22\lambda _2^26\lambda _2=0`$ and have $`\lambda _2=3(\sqrt{21}4)/5`$ and $`\lambda _2=3`$, respectively.

$``$The third, $`O_2`$, is associated with the lower branch of the same equation and has $`\lambda _2=3(\sqrt{21}+4)/5`$.

The other five solutions, denoted by $`X_k`$ ($`k=1,\mathrm{},5`$) do possess nodes.

(a) Two two-component solutions $`X_{\pm 1}`$, again identical to those found in the pure $`{}_{}{}^{3}P_{2}^{}`$ pairing problem, with $`\lambda _1=0`$ and $`\lambda _2=\pm 1`$.

(b) Three three-component solutions, $`X_2`$, $`X_3`$, and $`X_4`$, associated with the parabola $`\lambda _2=\lambda _1^2/2+1`$ and having $`\lambda _2=132\sqrt{35}`$, $`\lambda _2=3`$, and $`\lambda _2=13+2\sqrt{35}`$, respectively.

These general features of the spectrum of solutions of the $`{}_{}{}^{3}P_{2}^{}`$$`{}_{}{}^{3}F_{2}^{}`$ problem are expected to persist even if $`|\eta |`$ is not so small.

7. THE BATTLE BETWEEN PHASES

To complete the phase diagram of superfluid neutron matter, we need the gap values $`\mathrm{\Delta }_F`$ for the various phases. These are found through the third basic equation (15c), which involves $`J_5`$. To determine which phase wins the competition at a given temperature $`T`$, we compare the free-energy shifts

$$F_s=_0^g\mathrm{\Delta }_F^2(g^{})\frac{dg^{}}{(g^{})^2}$$
$`(25)`$
due to pairing in the corresponding superfluid states, where $`g`$ is the relevant coupling constant.

At low $`T`$, the only viable contestants are the solutions with nodeless order parameters, since the other solutions lie too high in energy. The separation between the two groups of states simply cannot be bridged if the value of $`|\eta |`$ remains rather small.

As documented in the last section, the parametric degeneracy of the $`{}_{}{}^{3}P_{2}^{}`$ pairing problem in the $`\lambda _1\lambda _2`$ plane, embodied in the spectral relation (26), is completely eradicated when the $`\eta `$-coupling is switched on. On the other hand, the energetic degeneracy between the different superfluid phases is only partially lifted. Specifically, the spectrum of pairing energies decays into several groups of nearly degenerate states: The $`O_1`$ and $`O_2`$ phases form the lowest-energy group, followed by the phases $`O_{\pm 3}`$ along with the one-component phase with $`M=0`$, and so on.

Figure 3. Temperature dependence of the splitting between the phase $`O_1`$ (or $`O_2`$) and the phases $`O_{\pm 3}`$ (or the one-component phase with $`M=0`$), in terms of the difference $`\delta \mathrm{\Delta }_F^2(T)`$ between the respective $`\mathrm{\Delta }_F^2(T)`$ values, measured relative to $`|\delta \mathrm{\Delta }_F^2(T=0)|`$. To set the energy scales involved, the latter quantity is around 2% of $`\mathrm{\Delta }_F^2`$ for $`\eta =0.3`$, while the gap parameter $`\mathrm{\Delta }_F(T=0)`$ itself reaches a maximum value near 0.44 MeV at $`k_F2.1`$ fm-1 when the pairing interaction is given by the Argonne $`v_{18}`$ potential and free normal-state single-particle energies are employed .

Raising the temperature from $`T=0`$ toward the critical temperature $`T_c`$, the splitting between the two groups lowest in energy shrinks, until, at $`T0.7T_c`$, their roles are interchanged and transitions occur. This behavior is quantified in Fig. 3. We observe that the analogy with the A-B phase transition in superfluid 3He is imperfect. Whereas the order parameter of the A-phase solution has nodes while that of the B-phase does not, the transition in neutron matter takes place between phases with nodeless order parameters.

8. CONCLUSIONS

We close with some details and caveats and examine some prospects for observable implications of our findings.

First, there is the matter of unresolvable fine structure in the phase diagram. Introduction of the $`{}_{}{}^{3}P_{0}^{}`$ and $`{}_{}{}^{3}P_{2}^{}`$ pairing channels would entail a further lifting of energy degeneracies, but current ignorance of the in-medium interaction precludes reliable estimation of these effects.

Secondly, we point out that the region near the critical temperature $`T_c`$ may admit subtle phase transitions not uncovered by our analysis. In this regime, strong-coupling corrections  can no longer be neglected in deciding the contest between phases. Also, the character of the phase diagram may be influenced significantly by external magnetic fields, a fact not considered in our treatment. And further, the phenomenon of fermion condensation , occurring as a precursor to pion condensation \[12-14\], may provide an exotic source of phase transitions in superfluid neuron matter.

In summary, it has been established that the superfluid phase diagram of dense neutron matter can in principle exhibit several triplet superfluid phases (at least 13!). Transitions between the different phases are expected to occur as a young neutron star cools. Since the gap value changes in these transitions, their occurrence may ultimately be detected in the thermal history and/or the rotational dynamics of the star. Such transitions may produce a significant variation of the moment of inertia of the star, or alterations of the distribution of the angular momentum between the crust and the vortex system, resulting in a change of the star’s angular velocity.

The existing analysis of triplet pairing in neutron matter provides a firm foothold for the more ambitious project of a complete characterization of the repertoire of superfluid phases of liquid 3He, which is even more complex, promising a richness of phenomena not yet revealed by experiment or theory. Since $`{}_{}{}^{3}P_{2}^{}`$, $`{}_{}{}^{3}P_{1}^{}`$, and $`{}_{}{}^{3}P_{0}^{}`$ states all enter the picture, one will be faced with nine nonlinear equations for the relevant gap components (in contrast to the five encountered in the pure $`{}_{}{}^{3}P_{2}^{}`$ case). Nevertheless, the separation method will continue to provide an incisive tool of analysis.

Thus, the procedure applied here may in principle be extended to count all possible solutions of the set of nine BCS equations for the gap function of superfluid 3He. We may begin by severing the couplings between $`{}_{}{}^{3}P_{2}^{}`$, $`{}_{}{}^{3}P_{1}^{}`$, and $`{}_{}{}^{3}P_{0}^{}`$ components, which facilitates a reference analytic solution. With this solution at hand, we may restore the channel couplings and solve the resulting set of equations numerically step by step, under increase of the coupling parameters. Attention then turns to (i) verification of the structure of the A-phase predicted by ABM , (ii) evaluation of the difference between free energies of A and B phases vs. temperature and pressure, and the barrier between the two phases (important for understanding the huge delay of the phase transition between the supercooled A-phase and the B-phase \[17–19\]), (iii) investigation of strong-coupling corrections beyond the Ginzburg-Landau approximation, and (iv) examination of the phase diagram of liquid 3He in an external magnetic field.

ACKNOWLEDGMENTS

The research described herein was supported by the U. S. National Science Foundation under Grant No. PHY-9900713 (JWC and VAK), by the McDonnell Center for the Space Sciences (VAK), and by Grant No. 00-15-96590 from the Russian Foundation for Basic Research (VAK and MVZ). JWC is grateful to the US Army Research Office–Research Triangle Park for travel support through a grant to Southern Illinois University–Carbondale. We thank A. Sedrakian, I. I. Strakovsky, G. E. Volovik, and D. N. Voskresensky for stimulating and informative discussions.

REFERENCES

 D. D. Osheroff, Rev. Mod. Phys. 69, 667 (1997).

 T. Takatsuka and R. Tamagaki, Prog. Theor. Phys. Suppl. 112, 27 (1993).

 V. A. Khodel, V. V. Khodel, and J. W. Clark, Nucl. Phys. A598, 390 (1996).

 V. A. Khodel, V. V. Khodel, and J. W. Clark, Phys. Rev. Lett. 81, 3828 (1998).

 J. W. Clark, V. A. Khodel, and V. V. Khodel, in Condensed Matter Theories, Vol. 13, ed. J. da Providencia and F. B. Malik (Nova Science Publishers, Commack, NY, 1998).

 T. Lindenau, M. L. Ristig, and J. W. Clark, in Condensed Matter Theories, Vol. 14, ed. D. J. Ernst, I. Perakis, and S. Umar (Nova Science Publishers, Commack, NY, 1999), pp. 131-139.

 J. W. Clark, V. V. Khodel, and V. A. Khodel, in Condensed Matter Theories, Vol. 15, ed. G. S. Anagnostatos, R. F. Bishop, K. A. Gernoth, J. Ginis, and A. Theophilou (Nova Science Publishers, Huntington, NY, 2000), pp. 1-12.

 V. V. Khodel, V. A. Khodel, and J. W. Clark, Nucl. Phys. A679, 827 (2001).

 D. Vollhardt and P. Wölfle, The Superfluid Phases of Helium 3 (Taylor & Francis, London, 1990).

 J. A. Sauls and J. W. Serene, Phys. Rev. D 17, 1524 (1978).

 V. A. Khodel and V. R. Shaginyan, JETP Lett. 51, 553 (1990).

 D. N. Voskresensky, V. A. Khodel, M. V. Zverev, and J. W. Clark, Ap. J. 533, L127 (2000).

 J. W. Clark, V. A. Khodel, and M. V. Zverev (2001) Physics of Atomic Nuclei (Yad. Fiz.) 64, 619-626 (2001).

 J. W. Clark, V. A. Khodel, and M. V. Zverev, in Condensed Matter Theories, Vol. 16, ed. S. Hernandez and J. W. Clark (Nova Science Publishers, Huntington, NY, 2001), pp. 263-273.

 P. W. Anderson and P. Morel, Phys. Rev. 123, 1911 (1961).

 P. W. Anderson and W. F. Brinkman, Phys. Rev. Lett. 30, 1108 (1973).

 A. J. Leggett, J. Low Temp. Phys. 87, 571 (1992).

 D. D. Osheroff and M. C. Cross, Phys. Rev. Lett. 38, 905 (1977).

 D. S. Buchanan, G. W. Swift, and J. C. Wheatley, Phys. Rev. Lett. 57, 341 (1986).