text:
# What ignites optical jets?To appear in the proceedings of the workshop “Relativistic jets in the Chandra and XMM Era” (New AR) This work was supported by the U.S. Department of Energy under contract No. DE-AC02-76CH03000.

## 1 Introduction

In the standard model for Active Galactic Nuclei (AGN), jets transfer mass, energy (kinetic and electromagnetic) and momentum from the central object into the surrounding medium. Liu and Zhang (2002) recently compiled a list of radio jets from the literature. Adopting the definition from Bridle and Perley (1984), they list 925 jets associated with 661 radio sources. Most of the jets have been identified by their radio synchrotron emission. Only a few of them also show optical synchrotron emission. Based on the most recent list of optical jets by Scarpa and Urry (2002) and the list of X-ray jets by Harris and Krawczynski (2002), I have compiled a list of 26 jets from 26 sources which have optical emission associated with them111The full list with references is available online at http://home.fnal.gov/~jester/optjets and will be contained in a forthcoming publication. The following optical jets are not in Scarpa and Urry (2002): 3C 31, PKS 0637$``$752, B2 0755+37, PKS 1136$``$135, PKS 1150+497, 3C 279, 3C 293, PKS 1354+195, 3C 303, B2 1553+254, B2 1658+30A, 3C 380..

Only a few optical jets have been studied in detail, most notably M87 (e.g., Meisenheimer et al., 1996; Perlman et al., 1999, 2001) and, in fewer bandpasses, 3C 273 (e.g., Röser and Meisenheimer, 1991; Jester, 2001, and references therein), but they have hardly been considered as a class of objects. The small number of known optical jets suggests the conclusion that they are very rare objects. This would, in fact, be expected from synchrotron physics: the magnetic fields expected in jets are of order of a few nanoTesla. This means that the loss timescales for electrons emitting optical synchrotron radiation are of order 1000 years or shorter, while typical arcsecond-scale jets have lengths from a few up to a few hundred kiloparsec. Thus, it seems natural that the frequency $`\nu _\mathrm{c}`$ at which a jet’s power-law spectrum cuts off lies below the optical wavelength region in most cases.

However, in most jets where optical emission has been observed, it extends over scales much larger than the lifetimes of “optical” electrons would permit. Scarpa and Urry (2002) showed that *on average*, optical jets are sufficiently relativistic to remove the lifetime discrepancy by time dilation, without the need for reacceleration of particles within the jets (before this, Sparks et al., 1995, had already shown that all optical jets known then are probably beamed). Still, at least in 3C 273’s jet in-situ acceleration is necessary to explain the observed optical synchrotron emission, in particular because very strong beaming increases the energy loss by inverse Compton scattering of microwave background photons (Jester et al., 2001). Moreover, the X-rays from many jets are explained as synchrotron emission as well (Harris and Krawczynski, 2002) – with much shorter loss timescales than in the optical, and correspondingly a more severe lifetime discrepancy which is harder to remedy by relativistic effects.

It appears thus that the number of optical and X-ray jets is small because some special physical conditions are necessary either to switch on the extended acceleration mechanism (whose nature remains elusive), or to accelerate the corresponding electrons near the AGN’s core and allow their escape to kiloparsec scales. However, in direct contradiction to this, Capetti et al. (2000) noted that the detection rate of optical jets in the B2 radio survey implies they are “not particularly rare”.

Here I make an attempt to find out what the required physical conditions are, and indeed whether they are special. Like Sparks et al. (1995), I find that beaming plays a role, and that more work is necessary before we can get beyond the observer-dependent beaming to the intrinsic physical conditions.222The cosmology is that used by Liu and Zhang (2002): $`H_0=100`$ km/s/Mpc, $`q_0=0.5`$.

## 2 A case study: 3C 273

In this place, I only briefly mention an idea which could help to explain the onset of bright high-frequency emission in the jet in 3C 273, where beaming alone cannot explain the observed optical emission. This jet has a divided appearance, with a continuous radio jet extending out to 21$`.^{\prime \prime }`$4 projected distance from the core, but with bright optical and X-ray emission only beyond 12′′ from the core. In Jester & Krause , we suggest that the jet may be lit up by making a transition from being overdense to underdense, with a strong shock forming near the transition region. We show a 2.5D hydrodynamical simulation supporting this suggestion. Full details will be contained in our forthcoming publication.

## 3 Optical jets: a beaming selection effect?

Synchrotron radio sources are usually assumed to have approximately power-law spectra, with a steep high-energy cutoff at some frequency $`\nu _\mathrm{c}`$ above which there is no significant emission. An objective definition of an “optical” jet would require a cutoff frequency beyond the near-infrared. A statistically complete survey for optical jets would include a surface brightness limit. Such a survey does not exist. I have therefore divided the radio jet list by Liu and Zhang (2002) into subsets: the 26 jets with detected optical emission, and the remainder as “radio-only” jets. This is in effect a qualitative definition of an “optical” jet.

As a third class, I consider X-ray jets (this class consists entirely of objects already included in the previous two classes). Seven radio jets show X-ray emission, but no optical, and thirteen jets have both optical and X-ray emission. I am not considering optical or X-ray emission from hot spots, only from the jet beams themselves.

I first consider whether the presence of optical emission from the jet is related to the core’s radio power. Figure 1 shows that at fixed redshift, nearly all optical jets are associated with radio sources which are more powerful than the median at that redshift. The same is true for X-ray jets. Conversely, among the sources at a fixed radio power, only the lower-redshift sources have optical or X-ray emission. Thus, there is no intrinsic correlation between core power and the presence of high-frequency emission.

The histograms in Figure 1 show that optical and X-ray jets are more likely to be 1-sided than an average radio jet, although this arises in part from the correlation between sidedness and jet power (Liu and Zhang, 2002). Jet sidedness is of course an indicator for relativistic beaming.

Another beaming indicator is the core prominence ratio $`R`$, the ratio of core power (affected strongly by beaming) to extended power (not affected strongly) at a fixed rest-frame frequency (e.g., Laing et al., 1999). Figure 2 shows the core prominence ratio of the three classes of jets as function of extended power. I have computed $`R`$ from the observed-frame values for core power at 5 GHz and total power at 1.4 GHz listed by Liu and Zhang (2002), assuming a flat core spectrum and $`f_\nu \nu ^1`$ for the extended component.

Scarpa and Urry (2002) have argued that the emission from optical jets is indeed on average beamed up. But Figures 1 and 2 suggest a stronger statement: First, all sources with detected radio jets are beamed. Secondly and more importantly, optical and X-ray jets are *more strongly* beamed than “radio-only” jets. By comparison of the jet points with the median relation for all radio sources (dashed line in Fig. 2), sources with radio jets in turn seem more strongly beamed than the radio source population as a whole. However, in interpreting this result, one should remember the inhomogeneity of the jet data set and the relatively crude spectral corrections.

Since the X-ray emission of some jets is explained as beamed inverse Compton scattering of microwave background photons (see the overview by Harris and Krawczynski, 2002), it is expected to find strong beaming for some of the X-ray jets. However, optical emission is usually shown or assumed to be synchrotron emission (with the possible exception of the ultraviolet excess observed from 3C 273’s jet; Jester et al., 2002), without reference to beaming. While there may be a physical connection between the jet’s bulk Lorentz factor and the intrinsic optical brightness (because more kinetic energy is available to be converted into relativistic particles, for example), I propose here that the sample of optical jets detected so far has overall no special physical significance, but they have only been detected because their intrinsic surface brightness has been beamed up sufficiently.

Just before the time of this writing, a paper appeared as preprint which shows that exactly this surface brightness selection effect can account for the optical non-detection of all but two jets in HST imaging of 57 sources from the B2 sample (Parma et al., 2003). If this holds for all other searches, it implies that optical jets are much more common in all extragalactic radio sources than has been appreciated so far, and they have not been detected in significant numbers simply because they have not been looked for with sufficient surface brightness sensitivity.

## 4 Conclusions: Finding more optical jets

The present comparison of the radio properties of sources with detected optical or X-ray jets to those with “radio-only” jets has not brought to light any connection between intrinsic physical properties and the presence of optical or X-ray emission. Instead, it is found that optical jets are more strongly beamed on average than other radio jets are. I therefore suggest that the entire sample of presently known optical jets is determined mainly by surface brightness selection effects. Like Parma et al. (2003) have done for the B2 sample, this suggestion has to be tested by comparing the optical surface brightness expected from an extrapolation of radio maps with a suitable power law ($`f_\nu \nu ^{0.7}`$ represents most jet spectra well) to the surface brightness limits of available optical observations. For those jets without suitably deep observations, the most efficient search method will employ observations in the near-infrared. In this way, jets with cutoffs just below the optical are still detectable. Positive detections will be followed up by observations in a second bandpass to determine the spectral index, and hence the detectability at higher frequencies. Only with a firm sample of detections *and* non-detections can we analyse the dependence of the jet’s synchrotron cutoff frequency $`\nu _\mathrm{c}`$ on other observables and thus address the question of what ignites optical jets.

I am grateful to Martin Hardcastle and Robert Laing for valuable comments.


format: text

text:
# Generalized Fibonacci numbers and dimer statistics

## Abstract

We establish new product identities involving the $`q`$-analogue of the Fibonacci numbers. We show that the identities lead to alternate expressions of generating functions for close-packed dimers on non-orientable surfaces.

Department of Physics, Northeastern University, Boston, Massachusetts 02115

The recursion relation

$$_{n+2}=_n+_{n+1}$$
(1)
with the initial values $`_0=_1=1`$ produces the well-known Fibonacci (1170-1250) numbers $`1,1,2,3,5,8,13,21,34,55,\mathrm{}`$ . The Fibonacci sequence can be generalized in a number of ways . One generalization is the $`q`$-analogue of the sequence defined by the recursion relation

$$_{n+2}(q)=_n(q)+q_{n+1}(q),$$
(2)
and the initial condition $`_0(q)=1`$ and $`_1(q)=q`$. Here $`q`$ is any real or complex number. This yields the sequence $`1,q,q^2+1,q^3+2q,q^4+3q^2+1,\mathrm{}`$ .

It can be readily established that the generating function of $`_n(q)`$ is

$$\frac{1}{1qss^2}=\underset{n=0}{\overset{\mathrm{}}{}}_n(q)s^n.$$
(3)
Writing $`q=xx^1`$ and partial fractioning the left-hand side of (3), one obtains the explicit expression

$$_n(q)=\frac{x^{n+1}+(1)^nx^{(n+1)}}{x+x^1},n=0,1,2,\mathrm{}.$$
(4)
We have following identities which we state as a theorem.

Theorem: (A) For $`M=`$ integers and $`\theta _m=m\pi /(2M+1)`$, we have

$$_M(q)\pm i_{M1}(q)=\underset{m=1}{\overset{M}{}}\left(q2i(1)^m\mathrm{cos}\theta _m\right),M=1,2,\mathrm{}.$$
(5)
(B) For $`M=`$ integers we have

$$_M(q)=\underset{m=1}{\overset{M}{}}\left(q+2i\mathrm{cos}\frac{m\pi }{M+1}\right),M=1,2,\mathrm{}.$$
(6)
(C) For $`N=`$ integers and $`\varphi _n=(2n1)\pi /2N`$, we have

$`_N(q)+_{N2}(q)`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{N}{}}}\left(q\pm 2i\mathrm{cos}\varphi _n\right)N=2,4,\mathrm{}`$ (7)
$`=`$ $`\pm 2i+{\displaystyle \underset{n=1}{\overset{N}{}}}\left(q\pm 2i(1)^n\mathrm{sin}\varphi _n\right),N=1,3,\mathrm{}.`$ (8)
(D) For $`M=`$ integers, $`N`$ = odd, $`q_n=2\tau ^1\mathrm{sin}\varphi _n`$, and $`p_m=2\tau \mathrm{cos}\theta _m`$, we have

$$\tau ^{MN}\underset{n=1}{\overset{N}{}}\left(_M(q_n)i(1)^{n+M}_{M1}(q_n)\right)=i^{M^2}\underset{m=1}{\overset{M}{}}\left(_N(p_m)+_{N2}(p_m)+2i(1)^{m+M}\right),$$
(9)
These identities can be established by determining the zeroes of the expressions. To establish $`(A)`$, for example, using (4) it is straightforward to deduce the identity

$$_M(q)\pm i_{M1}(q)=x^M\left(\frac{x^{2M+1}(\pm i)^{2M+1}}{xi}\right).$$
(10)
The quantity inside the parentheses in (10) is a polynomial of degree $`2M`$ in $`x`$ with the highest term $`x^{2M}`$. It has $`2M`$ zeroes at

$`x_m=\pm ie^{i\theta _{2m}},m=1,2,\mathrm{},2M.`$
Using the fact that $`x_{2M}=x_1^{},x_{2M1}=x_2^{},\mathrm{},x_{M+1}=x_M^{}`$, we obtain

$`_M(q)\pm i_{M1}(q)`$ $`=`$ $`x^M{\displaystyle \underset{m=1}{\overset{2M}{}}}(xx_m)`$ (11)
$`=`$ $`x^M{\displaystyle \underset{m=1}{\overset{M}{}}}(xx_m)(x+x_m^{})`$
$`=`$ $`{\displaystyle \underset{m=1}{\overset{M}{}}}\left(q2i\mathrm{cos}\theta _{2m}\right).`$
This leads to $`(A)`$ after we make use of the identities $`\mathrm{cos}\theta _{2M}=\mathrm{cos}\theta _1,\mathrm{cos}\theta _{2M2}=(1)^3\mathrm{cos}\theta _3,\mathrm{}`$. The identity $`(B)`$ is proved in a similar fashion.

To establish $`(C)`$ we note that for all $`N`$ we have

$`_N(q)+_{N2}(q)=x^N+(x)^N.`$
For $`N=`$ even, the polynomial $`x^N(x^{2N}+1)`$ has $`2N`$ zeroes at

$$x_n=ie^{i\varphi _n},x_n=ie^{i\varphi _n},n=1,2,\mathrm{},N.$$
(12)
This leads to (7) by factorizing the polynomial as in the proof of $`(A)`$.

For $`N=`$ odd, we have instead

$`_N(q)+_{N2}(q)2i`$ $`=`$ $`x^Nx^N2i`$ (13)
$`=`$ $`x^N\left(x^Ni\right)^2`$
$`=`$ $`x^N{\displaystyle \underset{n=1}{\overset{N}{}}}(xx_n)^2`$
where $`x_n=(1)^ne^{i\varphi _n}.`$ The identity (8) now follows the replacement of $`x_1=x_N^{},x_2=x_{N1}^{},\mathrm{}`$ in one of the $`xx_n`$ factors.

To establish $`(D)`$, write $`q=q_n=2\tau ^1\mathrm{sin}\varphi _n`$ in $`(A)`$ and $`q=p_m=2\tau \mathrm{cos}\theta _m`$ in $`(C)`$ to obtain, respectively,

$`_M(q_n)i(1)^{n+M}_{M1}(q_n)`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{M}{}}}\left(2\tau ^1\mathrm{sin}\varphi _n+2i(1)^{m+n+M}\mathrm{cos}\theta _m\right)`$
$`_N(p_m)+_{N2}(p_m)+2i(1)^{m+M}`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{N}{}}}\left(2\tau \mathrm{cos}\theta _m2i(1)^{n+m+M}\mathrm{sin}\varphi _n\right).`$
The identity $`(D)`$ now follows by combining these two expressions and the fact that, for $`N=`$ odd integers,

$`{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \underset{n=1}{\overset{N}{}}}\left(i(1)^{m+M+n+1}\right)=i^{M^2}.`$
The identities $`(A)`$ \- $`(D)`$ can be used to deduce alternate expressions for dimer generating functions on a Möbius strip. Let $`Z_{,𝒩}^{\mathrm{Mob}}(z_h,z_v)`$ be the generating function on an $`\times 𝒩`$ Möbius strip of width $`𝒩`$ with horizontal dimer weights $`z_h`$ and vertical dimer weights $`z_v`$. It can be shown that, by using $`(A)`$ \- $`(D)`$ and the equivalences of the sets $`\{\mathrm{sin}\varphi _{2n}\}=\{\mathrm{cos}\varphi _n\}`$ for $`N=`$ even and $`\{\mathrm{sin}\varphi _{2n}\}=\{(1)^n\mathrm{sin}\varphi _n\}`$ with $`n=1,2,\mathrm{},N`$, one establishes the following equivalent expressions:

For $`N=`$ even and $`\tau =z_v/z_h`$,

$`Z_{2M,N}^{\mathrm{Mob}}(z_h,z_v)`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \underset{n=1}{\overset{N/2}{}}}\left(4z_h^2\mathrm{sin}^2\varphi _{2n}+4z_v^2\mathrm{cos}^2\theta _m\right)`$ (14)
$`=`$ $`z_h^{MN}{\displaystyle \underset{m=1}{\overset{M}{}}}\left(_N(p_m)+_{N2}(p_m)\right)`$ (15)
$`=`$ $`z_v^{MN}{\displaystyle \underset{n=1}{\overset{N/2}{}}}_{2M}(q_n),N=\mathrm{even}`$ (16)
$`Z_{2M1,N}^{\mathrm{Mob}}(z_h,z_v)`$ $`=`$ $`2z_h^{N/2}{\displaystyle \underset{m=1}{\overset{M1}{}}}{\displaystyle \underset{n=1}{\overset{N/2}{}}}\left(4z_h^2\mathrm{sin}^2\varphi _{2n}+4z_v^2\mathrm{cos}^2(m\pi /2M)\right)`$ (17)
$`=`$ $`2z_h^{(2M1)N/2}{\displaystyle \underset{m=1}{\overset{M1}{}}}[_N\left(2\tau \mathrm{cos}(m\pi /2M)\right)+_{N2}\left((2\tau \mathrm{cos}(m\pi /2M))\right]`$ (18)
$`=`$ $`2z_h^{N/2}z_v^{(M1)N}{\displaystyle \underset{n=1}{\overset{N/2}{}}}\left(q_n^1_{2M1}(q_n)\right),N=\mathrm{even}.`$ (19)
For $`N=`$ odd and $`\tau =z_v/z_h`$,

$`Z_{2M,N}^{\mathrm{Mob}}(z_h,z_v)`$ $`=`$ $`\mathrm{Re}\left[(1i){\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \underset{n=1}{\overset{N}{}}}\left(2i(1)^{M+m+1}z_h\mathrm{sin}\varphi _{2n}+2z_v\mathrm{cos}\theta _m\right)\right]`$ (20)
$`=`$ $`z_h^{MN}\mathrm{Re}\left[(1i){\displaystyle \underset{m=1}{\overset{M}{}}}\left(_N(p_m)+_{N2}(p_m)+2i(1)^{m+M}\right)\right]`$ (21)
$`=`$ $`z_v^{MN}\mathrm{Re}\left[(1i)i^{M^2}{\displaystyle \underset{n=1}{\overset{N}{}}}\left(_M(q_n)i(1)^{n+M}_{M1}(q_n)\right)\right],N=\mathrm{odd}.`$ (22)
Expressions (14), (17) and (20) are those given by Lu and Wu , expressions (15), (18) and (21) are given by Tesler , and expressions (16), (19) and (22) are new.

Work has been supported in part by NSF grant DMR-9980440.


text:
# 1 Introduction

## 1 Introduction

The Skyrme model is a non-linear field theory in which nuclei are modelled by quantum states of classical topological solitons . In this talk I will use the term topological soliton for minimal energy solutions in a field theory which have an associated integer degree or topological charge. The degree is conserved for topological reasons and ensures the stability of the soliton. In the Skyrme model static configurations are maps

$`U:^3SU(2).`$ (1.1)
Points in $`^3`$ will be denoted by $`x=(x_1,x_2,x_3)`$ and I write $`r`$ for the Euclidean length $`|x|=\sqrt{x_1^2+x_2^2+x_3^2}`$. Physically, $`U`$ combines the pion fields $`\pi _1,\pi _2,\pi _3`$ and the $`\sigma `$-field

$`U(x)=\sigma (x)+i\pi _a(x)\tau _a,`$ (1.2)
where $`\tau _a`$, $`a=1,2,3`$, are the Pauli matrices and the constraint $`\sigma ^2+\pi _1^2+\pi _2^2+\pi _3^2=1`$ is imposed. Here and in the following, summation over repeated indices is always implied. The Skyrme energy functional is

$`E[U]={\displaystyle d^3x\left(\frac{1}{2}\mathrm{tr}(L_iL_i)+\frac{1}{16}\mathrm{tr}([L_j,L_i][L_j,L_i])\right)},`$ (1.3)
where

$`L_i=U^{}_iU`$ (1.4)
and $`_i=/x_i`$ with $`i=1,2,3`$. The Euler-Lagrange equation can be expressed in terms of the modified currents

$`\stackrel{~}{L}_i=L_i{\displaystyle \frac{1}{4}}[L_j,[L_j,L_i]]`$ (1.5)
and reads

$`_i\stackrel{~}{L}_i=0.`$ (1.6)
It was already noted by Skyrme that the finite-energy requirement means that the Skyrme fields have to tend to a constant at spatial infinity, thus becoming effectively maps from $`S^3`$ to $`SU(2)S^3`$ with an associated integer degree. Skyrme identified the degree physically with the baryon number. The first rigorous proof that for finite-energy Skyrme configuration the degree

$`\text{deg}[U]={\displaystyle \frac{1}{24\pi ^2}}{\displaystyle d^3xϵ_{ijk}\mathrm{tr}\left(L_iL_jL_k\right)}`$ (1.7)
is an integer was only given fairly recently in . As a result the configuration space

$`𝒞=\{U:^3SU(2)|E[U]<\mathrm{}\}`$ (1.8)
is partitioned into sectors $`𝒞_k`$ consisting of all finite energy configurations of degree $`k`$. Faddeev’s bound

$`E[U]>12\pi ^2|k|`$ (1.9)
(the strict inequality follows from ) implies the existence of infima

$`I_k=\text{inf}\{E[U]|U𝒞_k\}.`$ (1.10)
The central question of this talk is whether the infima are attained, i.e. whether minima exist in all sectors of the Skyrme model. In the following I call minimal energy solutions of non-vanishing degree Skyrme solitons.

To illustrate quite how little is known rigorously about the existence of Skyrme solitons note that the hedgehog ansatz

$`U_H(x)=\mathrm{exp}(if(r)\widehat{x}_a\tau _a)`$ (1.11)
with the boundary condition $`f(0)=\pi `$ and $`f(\mathrm{})=0`$ leads to a minimisation problem for the profile function $`f`$, which is known to have a unique solution , called the Skyrmion. It is also known  that the Skyrme energy functional has a minimum in the sector $`𝒞_1`$. However, it is still not known whether the minimum is of the spherically symmetric hedgehog form (1.11).

Physically, we may think of Skyrme solitons of degree $`k>1`$ as bound states of Skyrmions. One therefore expects the questions of whether such bound states exist to be related to the existence of attractive forces in the Skyrme model. This expectation is borne out by the analytical work of Esteban in . Esteban showed that for a suitable class functions

$`I_kI_l+I_{kl}`$ (1.12)
for all $`k,l`$. Physically this result is only sufficient to ensure threshold bound states \- it would be satisfied if for example $`I_k=C|k|`$ for some constant $`C`$. Esteban went on to show that a minimum exist in $`𝒞_k`$ provided one assumes the strict inequality

$`I_k<I_l+I_{kl}`$ (1.13)
for all $`k\{0,\pm 1\}`$ and $`l\{0,k\}`$ in the range $`|l|+|kl|<\sqrt{2}|k|`$. This result will be referred to as Esteban’s theorem in the following discussion. In  it was shown that the result still holds if one widens the class of allowed functions, but the inequality (1.13) remains a necessary assumption in the proof. In the cases where the infima $`I_l`$ and $`I_{kl}`$ are attained by Skyrme solitons, the inequality (1.13) is equivalent to the existence of attractive forces between those Skyrme solitons. In  the existence of attractive forces is studied in detail for Skyrme solitons obeying certain regularity assumptions. In the second part of this talk I review the arguments and results of that paper and in the third part I discuss their relation to Esteban’s result.

Before plunging into the technical analysis I should point out that much is known numerically about Skyrme solitons. Numerical searches aided by analytical ansätze suggest the existence of Skyrme solitons for $`1k22`$ . However, it is also worth stressing that the existence of classical minima does not ensure the existence of a quantum bound state. This is relevant in the physical application of the Skyrme model because nuclei do not exist for arbitrarily large baryon number. The study of the deuteron as a quantum state of the $`k=2`$ toroidal Skyrme soliton in  furthermore shows that when such a quantum bound state exists the quantum mechanical matter distribution may be quite different from that of the underlying classical soliton.

## 2 The interaction energy of Skyrme solitons

### 2.1 Asymptotics of a Skyrme soliton

The goal of this subsection is to study the behaviour of Skyrme solitons near spatial infinity. The key step is to think of the 2-sphere at infinity as a boundary of $`^3`$ and to show that the Skyrme equation is regular there. It then follows from a unique continuation argument that non-trivial solutions of the Skyrme equation have a non-trivial expansion near infinity. In other words, Skyrme solitons necessarily have a non-trivial large $`r`$ expansion in powers of $`1/r`$, possibly combined with $`\mathrm{ln}r`$.

To explain the basic ideas, consider a finite point $`x_0`$ and assume (if necessary redefining $`U(x)U(x_0)^1U(x)`$ ) that $`U(x_0)=1`$. Then we have the expansion

$`U(x)=1+u(x)`$ (2.1)
in a neighbourhood of $`x_0`$, with the $`2\times 2`$ complex matrix $`u`$ satisfying the algebraic constraints

$`u+u^{}+uu^{}=0,\mathrm{tr}(u)+\text{det}(u)=0.`$ (2.2)
After lengthy algebra, the Skyrme equation can be rewritten

$`P(u,u,^2u)=Q(u,u)+F(u,u),`$ (2.3)
where $`P`$ is linear in the second derivatives $`^2u`$ and contains the Laplace operator $`\mathrm{\Delta }u`$. The functions $`Q`$ and $`F`$ crucially do not depend on $`^2u`$ and are, respectively, quadratic and of degree four in $`u`$. The point of writing the Skyrme equation in this way is that the operator

$`fP(u,u,^2f)`$ (2.4)
is linear and elliptic near $`x_0`$. If we assume that both $`u`$ and the currents $`L_i`$ are Hölder continuous, we deduce by elliptic regularity that a solution of (2.3) is twice differentiable. Bootstrapping further, one concludes that any solution is in fact smooth (i.e. $`C^{\mathrm{}}`$). The key idea in this procedure is to re-interpret the non-linear Skyrme equation as an elliptic linear equation with coefficients depending on $`u`$ and $`u`$.

In order to apply the same idea at spatial infinity, we treat the 2-sphere at spatial infinity as a boundary of $`^3`$ and introduce coordinates $`(s,\theta ,\phi )`$, where $`s=1/r`$ and $`\theta `$ and $`\phi `$ are the usual spherical coordinates. Then defining

$`D_i=r_i={\displaystyle \frac{1}{s}}_i`$ (2.5)
the Euclidean Laplacian takes the form

$`\mathrm{\Delta }=s^2\mathrm{\Delta }_b`$ (2.6)
with

$`\mathrm{\Delta }_b=s^2_s^2+\mathrm{\Delta }_\omega ,`$ (2.7)
where $`\mathrm{\Delta }_\omega `$ is the Laplacian of the unit 2-sphere. Now we write $`U=1+u`$ for large $`r`$ (that is, for small positive $`s`$) and obtain a ’$`b`$’-version of the Skyrme equation

$`P_b(u,Du,D^2u)=Q_b(u,Du)+s^2F_b(u,Du).`$ (2.8)
This equation can now be analysed using tools from the theory of so-called $`b`$-differential operators, see . In particular one can use the boundary conditions that $`U`$ tends to one at $`\mathrm{}`$ and that the currents $`L_i`$ vanish there to solve (2.8) iteratively, thus producing an expansion of the solution in powers of $`s=1/r`$, with possible factors $`\mathrm{ln}r`$ in higher order terms. At this stage we have an asymptotic expansion of Skyrme solitons, but do not yet know whether there is a leading non-zero term in the expansion. This can be established with the help of a unique continuation theorem , suitably adapted for our purposes.

The upshot of this chain of arguments, described in detail in , is that every Skyrme soliton has a leading Lie-algebra valued multipole field

$`u_M(x)=i\tau _a{\displaystyle \underset{m=M}{\overset{M}{}}}{\displaystyle \frac{4\pi }{2M+1}}Q_{Mm}^a{\displaystyle \frac{Y_{Mm}(\theta ,\phi )}{r^{M+1}}},`$ (2.9)
where $`Y_{Mm}`$ are the usual spherical harmonics on $`S^2`$. In accordance with the usual nomenclature we refer to the multipole field (2.9) as a $`2^M`$-pole and call $`M`$ the order of the multipole. The leading multipole moments $`Q_{Mm}^a`$ are independent of the location of the Skyrme soliton, and are acted on naturally by rotations and iso-rotations. They are crucial for the calculations in the following section. Note that it was already shown in  that the leading multipole cannot be a monopole. The Skyrmion field (1.11) is known to have an iso-triplet of dipoles as leading multipoles. For other Skyrme solitons, too, the leading multipoles have been investigated to some extent. The leading multipole of the largest known order is an octupole (i.e. multipole order 3) which is believed to arise in the charge seven icosahedral Skyrme soliton

### 2.2 Interaction energy of two scalar multipoles

As a technical preparation for the computation of the interaction energy of two Skyrme solitons we derive a formula for the interaction energy of two scalar multipoles. In our description of multipoles we use the conventions of  throughout. Consider the field of a $`2^M`$-pole centred at $`X_+=(0,0,R/2)`$, where $`R>0`$. In terms of spherical coordinates $`(\theta _+,\phi _+)`$ centred at $`X_+`$ it reads

$`f_M(x)={\displaystyle \frac{4\pi }{(2M+1)}}{\displaystyle \underset{m=M}{\overset{M}{}}}Q_{Mm}{\displaystyle \frac{Y_{Mm}(\theta _+,\phi _+)}{|xX_+|^{M+1}}}.`$ (2.10)
Similarly define $`X_{}=(0,0,R/2)`$ and let $`(\theta _{},\phi _{})`$ be spherical coordinates centred at $`X_{}`$. A $`2^N`$-pole field centred at $`X_{}=(0,0,R/2)`$ has the form

$`g_N(x)={\displaystyle \frac{4\pi }{(2N+1)}}{\displaystyle \underset{n=N}{\overset{N}{}}}\stackrel{~}{Q}_{Nn}{\displaystyle \frac{Y_{Nn}(\theta _{},\phi _{})}{|xX_{}|^{N+1}}}.`$ (2.11)
A natural measure for the interaction energy between harmonic function in the upper and those in the lower half plane is

$`V[f,g]={\displaystyle _{x_3=0}}𝑑x_1𝑑x_2(g_3ff_3g).`$ (2.12)
The computation of this interaction energy for the multipole fields $`f_M`$ and $`g_N`$ is not easy because the combined field of the two multipoles only has cylindrical symmetry about the $`x_3`$-axis and not the full rotational symmetry of each of the multipoles. As explained in , Fourier transformation in the $`x_1x_2`$-plane turns out to be an efficient tool for evaluating the integral in (2.12). Assuming without loss of generality that $`MN`$, the answer is

$`V[f_M,g_N]={\displaystyle \frac{(4\pi )^2}{R^{M+N+1}}}`$ $`{\displaystyle \frac{(M+N)!(1)^{N+M}}{\sqrt{(2M+1)(2N+1)}}}\times `$ (2.13)
$`{\displaystyle \underset{m=M}{\overset{M}{}}}{\displaystyle \frac{\overline{Q}_{Mm}\stackrel{~}{Q}_{Nm}}{\sqrt{(Mm)!(M+m)!}\sqrt{(Nm)!(N+m)!}}}.`$
For us, the most important feature of the formula (2.13) is that, for non-vanishing multipole moments, the interaction energy can always be made non-zero by a suitable rotation of the multipole of the highest order, in our case $`N`$ (if $`M=N`$ it does not matter which of the multipoles gets rotated). By definition the multipole components $`\stackrel{~}{Q}_{Nn}`$ are the components of an element $`\stackrel{~}{Q}`$ of the $`(2N+1)`$-dimensional irreducible unitary representation $`W_N`$ of $`SO(3)`$, with $`GSO(3)`$ acting via $`\stackrel{~}{Q}_{Nn}_{n^{}=N}^NU_{nn^{}}^N(G)\stackrel{~}{Q}_{Nn^{}}`$. Think of the pairing of the multipole moments in (2.13) as a linear form

$`F_Q:W_N,\stackrel{~}{Q}{\displaystyle \underset{m=M}{\overset{M}{}}}{\displaystyle \frac{\overline{Q}_{Mm}\stackrel{~}{Q}_{Nm}}{\sqrt{(Mm)!(M+m)!}\sqrt{(Nm)!(N+m)!}}}.`$ (2.14)
Then it follows from the irreducibility of $`W_N`$ that $`U^N(G)\stackrel{~}{Q}`$ cannot lie in the kernel of $`F_Q`$ for all $`G`$. Therefore $`F_Q(U^N(G)\stackrel{~}{Q})0`$ for some $`G`$. Hence the interaction energy (2.13) is non-vanishing after rotating $`\stackrel{~}{Q}`$ with that $`G`$.

### 2.3 Interaction energy of well-separated Skyrme solitons

Consider now two Skyrme solitons $`U^{(1)}`$ and $`U^{(2)}`$ of degrees $`k`$ and $`l`$. Since the total energies of both $`U^{(1)}`$ and $`U^{(2)}`$ are finite there must be balls $`B_1`$ and $`B_2`$ in $`^3`$ so that most of the energy of $`U^{(1)}`$ and $`U^{(2)}`$ is concentrated in, respectively, $`B_1`$ and $`B_2`$. Outside the balls $`B_1`$ and $`B_2`$ the asymptotic analysis of the previous section applies. Suppose that the leading multipole of $`U^{(1)}`$ is a $`2^M`$-pole and the leading multipole of $`U^{(2)}`$ is a $`2^N`$-pole. Denoting the radii of $`B_1`$ and $`B_2`$ by $`D_1`$ and $`D_2`$ we have

$`U^{(1)}(x)1+u_M(x)\text{for}xB_1`$ (2.15)
and

$`U^{(2)}(x)1+v_N(x)\text{for}xB_2,`$ (2.16)
where $`u_M`$ is of the form (2.9) and analogously

$`v_N(x)=i\tau _a{\displaystyle \underset{n=N}{\overset{N}{}}}{\displaystyle \frac{4\pi }{2N+1}}\stackrel{~}{Q}_{Nn}^a{\displaystyle \frac{Y_{Nn}(\theta ,\phi )}{r^{N+1}}}.`$ (2.17)
Using the translational invariance of the Skyrme energy functional we can assume without loss of generality that $`B_1`$ is centred at $`X_+=(0,0,R/2)`$ and that $`B_2`$ is centred at $`X_{}=(0,0,R/2)`$, where $`R`$ is so large that $`B_1`$ and $`B_2`$ do not overlap, i.e. $`R>D_1+D_2`$. The parameter $`R`$ will be interpreted as the separation of the Skyrme solitons. Then we define the following product configuration

$`U_R(x)=U^{(1)}(x)U^{(2)}(x).`$ (2.18)
This configuration has degree $`k+l`$ and finite energy, so that $`U_R𝒞_{k+l}`$. A lengthy calculation performed in  shows that

$`E[U_R]=E[U^{(1)}]+E[U^{(2)}]+\mathrm{\Delta }E+𝒪\left({\displaystyle \frac{1}{R^{2N+4}}}\right)+𝒪\left({\displaystyle \frac{1}{R^{2M+4}}}\right),`$ (2.19)
where

$`\mathrm{\Delta }E=2{\displaystyle \underset{a=1}{\overset{3}{}}}{\displaystyle _{x_3=0}}𝑑x_1𝑑x_2(u_M^a_3v_N^av_N^a_3u_M^a).`$ (2.20)
Now we note that $`\mathrm{\Delta }E`$ is the sum over iso-components of the the scalar interaction terms we studied above

$`\mathrm{\Delta }E_a=2{\displaystyle _{x_3=0}}𝑑x_1𝑑x_2(u_M^a_3v_N^av_N^a_3u_M^a)=2V[u_M^a,v_N^a].`$ (2.21)
Picking one of the iso-indices, say $`a=1`$, we can use iso-rotations to make sure that the first iso-components $`u_M^1`$ and $`v_N^1`$ are non-vanishing. The result of the previous subsection then implies that we can make the multipole interaction energy $`\mathrm{\Delta }E_1`$ non-zero by a spatial rotation of the Skyrme soliton with the higher multipole order.

Now consider the sum

$`\mathrm{\Delta }E=\mathrm{\Delta }E_1+\mathrm{\Delta }E_2+\mathrm{\Delta }E_3.`$ (2.22)
Following an idea in  we would like to show that we can always arrange for $`\mathrm{\Delta }E`$ to be negative by a suitable iso-rotation of one of the Skyrme solitons. We may assume that, possibly after re-labelling the pion fields,

$`\mathrm{\Delta }E_1\mathrm{\Delta }E_2\mathrm{\Delta }E_3.`$ (2.23)
If $`\mathrm{\Delta }E<0`$ we are done, so suppose that $`\mathrm{\Delta }E0`$. Since we know that not all $`\mathrm{\Delta }E_a`$ vanish we can conclude that $`\mathrm{\Delta }E_1>0`$. Now perform an iso-rotation of Skyrme soliton $`2`$ by $`180`$ degrees around the third iso-spin axis. This reverses the sign of $`v_N^1`$ and $`v_N^2`$ and hence of $`\mathrm{\Delta }E_1`$ and $`\mathrm{\Delta }E_2`$. The new value of $`\mathrm{\Delta }E`$ is

$`\mathrm{\Delta }E`$ $`=`$ $`\mathrm{\Delta }E_1\mathrm{\Delta }E_2+\mathrm{\Delta }E_3`$ (2.24)
$`=`$ $`\mathrm{\Delta }E_1(\mathrm{\Delta }E_2\mathrm{\Delta }E_3)<0,`$
since $`\mathrm{\Delta }E_1<0`$ and, with our ordering, $`(\mathrm{\Delta }E_2\mathrm{\Delta }E_3)0`$.

Thus, the contribution $`\mathrm{\Delta }E`$ to the interaction energy of two Skyrme solitons $`U^{(1)}`$ and $`U^{(2)}`$ can always be made negative by suitable relative rotations and iso-rotations of the Skyrme solitons. Returning to the expression (2.19) for the interaction energy and noting from (2.13) that $`\mathrm{\Delta }E`$ falls off like $`R^{(N+M+1)}`$, we conclude that $`\mathrm{\Delta }E`$ is the leading contribution for large separation $`R`$ provided the orders of the leading multipoles do not differ by more than 2, i.e. provided that $`|NM|2`$. In that case we can therefore always arrange for the interaction energy to be negative for sufficiently large separation. In symbols, for sufficiently large $`R`$ and the appropriate orientations we have

$`E[U_R]<E[U^{(1)}]+E[U^{(1)}].`$ (2.25)
## 3 The existence of Skyrme solitons

The results of the previous section are not sufficient to derive the existence minima in a general sector $`𝒞_k`$ of the Skyrme model. Nonetheless it is instructive to see how far one can get with the following two extra assumptions.

1. When minima exist they satisfy the technical assumptions made in . In particular the field $`U`$ and the currents $`L_i`$ are Hölder continuous.
2. The leading multipole of any Skyrme soliton is at most an octupole (order 3).

The second assumption seems particularly restrictive and deserves a comment. Since monopoles (order 0) cannot arise in Skyrme solitons and dipoles (order 1) are known to arise in some, the restriction to leading multipoles of order at most 3 means that any two Skyrme solitons we consider in the following have multipoles whose orders do not differ by more than 2. As mentioned earlier, none of the numerically studied Skyrme soliton violate assumption 2.

Recall from section 1 that Skyrme solitons have rigorously been shown to exist in the topological sectors $`𝒞_1`$ and $`𝒞_1`$. Although it is not clear whether the minimum is of the hedgehog form (1.11), we can use the numerically computed energy of the Skyrmion to give an upper bound on the infimum $`I_1`$. To avoid the discussion of numerical accuracy we use the energy of the instanton generated hedgehog field  $`\stackrel{~}{I}_1=1.24\times 12\pi ^2`$ to estimate

$`I_1<{\displaystyle \frac{5}{4}}\mathrm{\hspace{0.17em}12}\pi ^2.`$ (3.1)
Using Esteban’s weak inequality (1.12) we deduce

$`I_k<{\displaystyle \frac{5}{4}}\mathrm{\hspace{0.17em}12}\pi ^2|k|`$ (3.2)
for all $`k.`$ Combining this result with the energy bound (1.9) we also deduce that

$`I_k<{\displaystyle \underset{s}{}}I_{l_s}`$ (3.3)
for all $`k`$ and $`l_s`$, $`s`$ satisfying $`_sl_s=k`$ and $`_s|l_s|\frac{5}{4}|k|`$. Thus we can weaken the assumption in Esteban’s theorem about the existence of minima. The existence of a minimum in sector $`𝒞_k`$ follows if the strict inequality (1.13) is satisfied for all $`k`$ and all $`l\{k,0\}`$ satisfying

$`|l|+|kl|<{\displaystyle \frac{5}{4}}|k|.`$ (3.4)
Combining this result with our proof of attraction between Skyrme solitons and the assumptions made above we can prove the existence of minima in the sectors $`𝒞_k`$ for low values of $`k`$. We prove the claim by induction and first consider positive $`k`$. We know from  that the minimum exists for $`k=1`$. Suppose that minima exist in $`𝒞_l`$ with $`l=1,2,\mathrm{},k1`$. Then apply the product ansatz $`U_R`$ and the inequality (2.25) to pairs of Skyrme solitons $`U^{(1)}`$ and $`U^{(2)}`$ of degree $`l`$ and $`kl`$ for $`l=1,\mathrm{},k1`$. By assumption the energies of these Skyrme solitons are equal to the infima $`I_l`$ and $`I_{kl}`$. Then, by definition of the infimum $`I_k`$ (1.10) and from (2.25)

$`I_kE[U_R]<E[U^{(1)}]+E[U^{(1)}]=I_l+I_{(kl)}`$ (3.5)
for all $`l=1,\mathrm{},k1`$. As long as $`k8`$ the only integers $`l`$ satisfying the inequality (3.4) are in the range $`l=1,\mathrm{},k1`$, so that we have satisfied the condition of Esteban’s theorem and can conclude that the infimum in $`𝒞_k`$ is attained. The claim for negative $`k`$ follows by applying the energy preserving reflection map $`U(x)U^{}(x)`$.

For $`k>8`$ the step in the inductive proof fails because we then need to consider the inequality (1.13) for negative $`l`$, too. For example for $`k=9`$, the condition (3.4) is satisfied for $`l=1`$ and $`kl=10`$. However, since we do not know whether Skyrme solitons exist in the sector $`𝒞_{10}`$ we cannot establish the inequality $`I_9<I_1+I_{10}`$ with our method. Numerical evidence and physical intuition suggests that the inequality (1.13) should certainly hold when $`k`$ is positive and $`l`$ or $`(kl)`$ is negative. If it was violated in that range it would be energetically favourable for matter (solitons of positive degree) and antimatter (solitons of negative degree) to coexist in some sectors of the Skyrme model. However, while this seems unlikely, we cannot disprove the possibility with our methods. If one could rule it out by other methods, the inductive argument given above could be extended to prove the existence of minima in all sectors of the Skyrme model.

To end, I point out that it seems plausible that our results can also be used to give a more direct proof of the existence of minima in all sectors of the Skyrme model. The idea of such a proof is that the infimum in one of the sectors, say $`𝒞_k`$, can only fail to be attained if it is energetically favourable for configurations in that sector to split into well-separated Skyrme solitons. The infimum $`I_k`$ would then be realised by a virtual solution, made up of infinitely separated Skyrme solitons. The number of Skyrme solitons in this virtual solution is necessarily finite. For suppose the degrees of the Skyrme solitons are $`l_s`$, $`s`$ with $`_sl_s=k`$. By the inequality (3.3) we can only have the equality

$`I_k={\displaystyle \underset{s}{}}I_{l_s}`$ (3.6)
if $`_s|l_s|<\frac{5}{4}|k|`$, which means that all but a finite number of the $`l_s`$ must be zero. Suppose that $`l_s0`$ for $`s=1,2\mathrm{}S`$ and $`l_s=0`$ for $`s>S`$. Physically, one expects the following dichotomy to hold.

1. For given $`k`$ the infimum $`I_k`$ is attained by a Skyrme soliton of degree $`k`$.
2. The infimum $`I_k`$ is equal to a finite sum

$`I_k={\displaystyle \underset{s=1}{\overset{S}{}}}E_{l_s}`$ (3.7)
of energies $`E_{l_s}=I_{l_s}`$ of $`S`$ Skyrme solitons of degrees $`l_s\{0,k\}`$.

With our assumption on the leading multipoles our earlier proof of attractive forces in the Skyrme model implies that we can always lower the energy of a finite number of infinitely separated Skyrme solitons by chosing appropriate orientations for two of the Skyrme solitons and bringing them closer together. Thus the infimum $`I_k`$ cannot be of the form (3.7) and if the dichotomy is valid, option (a) must hold.

## 4 Conclusion

In this talk I reviewed the result, obtained jointly with Manton and Singer, that arbitrary Skyrme solitons attract for suitable relative orientation and iso-orientation, provided the orders of their leading multipole moments do not differ by more than two. I also indicated how this result, combined with the work of Esteban, can be used to shed light on the question of minima in the Skyrme model. However, the approaches sketched here also show how difficult it is to establish rigorous results on minima. In particular the assumption on the multipoles of Skyrme solitons severely restricts the generality of the sort of argument outlined above.

Acknowledgments I thank Nick Manton and Michael Singer for collaboration on the related paper  and discussions. I acknowledge financial support trough an EPSRC advanced research fellowship.


text:
# Laser cooling of a trapped two-component Fermi gas

## Abstract

We study the collective Raman cooling of a trapped two-component Fermi gas using quantum master equation in the festina lente regime, where the heating due to photon reabsorption can be neglected. The Monte Carlo simulations show, that 3D temperatures of the order of $`0.008`$ $`T_F`$ can be achieved. We analyze the heating related to background losses, and show that our laser-cooling scheme can maintain the temperature of the gas without significant additional losses.

The realization of Bose-Einstein condensation (BEC) BEC  and degenerated Fermi gases Fermions  has outbursted the interest in the physics of ultracold bosonic and fermionic gases. A Fermi gas with attractive interactions undergoes for temperatures $`T`$ below a critical one, $`T_c`$, a transition into the superfluid Bardeen-Cooper-Schrieffer (BCS) phase DeGennes . The accomplishment of BCS and its possible detection have been considered in detail (see e.g. Refs. BCS ). Unfortunately, evaporative cooling currently employed in experiments is based on collisional processes, and has not yet allowed to reach $`T<T_c`$, since $`T_c`$ is much smaller than the Fermi temperature $`T_F`$, and for $`T<T_F`$ the collisions are strongly suppressed due to Pauli blocking HollandF . Recently, however, it has been proposed that $`T_c`$ may be significantly increased by employing Feshbach resonances Holland , or optical lattices OptLatt . This opens the possibility to achieve BCS at a temperature where the collisions are still efficient enough.
In this Rapid Communication we study the laser cooling of a two-component Fermi gas in the festina lente (FL) regime Festina , where the spontaneous emission rate $`\gamma `$ is smaller than the trap frequency $`\omega `$. In this regime the heating due to photon reabsorption is prevented, as recently observed for tightly bound atoms in optical lattices Weiss . The laser cooling in the FL regime has been already predicted to work for bosons CoolBosons  and for polarized fermions fermicool . By using Monte Carlo simulations, we show here that laser cooling allows for bringing a two component Fermi system to $`TT_F`$.
The laser cooling of fermions toward $`TT_F`$ is obstacled by several problems. One of them is the inhibition of spontaneous emission, which results in the decrease of the cooling efficiency. In Ref. fermicool  we have predicted that the inhibition problem can be overcome by either dynamically adjusting the spontaneous emission rate in a Raman cooling process, or by employing specially designed anharmonic traps. In this paper we apply the former solution. Another obstacle is related to the inelastic losses that create holes deeply in the Fermi sea, producing a significant heating. We show that background collisions do not affect significantly the laser cooling. In fact, our cooling method can be employed to maintain a two-component Fermi gas at a fixed $`T`$, for a relatively long time in a trap. Finally, we derive an analytic formula that describes the heating of the trapped gas due to background collisions. The heating is shown to be smaller than in the homogeneous case.
We consider fermionic atoms with an accessible electronic three-level $`\mathrm{\Lambda }`$ scheme, containing states $`|g`$, $`|e`$ and $`|r`$. The ground state $`|g`$ is coupled via a Raman transition to the metastable state $`|e`$. The latter state is also coupled by an optical transition to the upper state $`|r`$, from which atoms rapidly decay into $`|g`$. The adiabatic elimination of $`|r`$, leads to an effective two-level system, characterized by tunable parameters: Rabi frequency $`\mathrm{\Omega }`$, and spontaneous emission rate $`\gamma `$. The value of $`\gamma `$ can be controlled by modifying the coupling from $`|e`$ to $`|r`$.

The atoms are placed in a dipole trap characterized by a Lamb-Dicke parameter $`\eta =2\pi a/\lambda `$, with $`a=\sqrt{\mathrm{}/2m\omega }`$ being the size of the ground state of the trap, and $`\lambda `$ the laser wavelength. A dipole trap was recently used for the all-optical production of a degenerate gas of two Li species Granade , and in current experiments in Mg Ernst . We consider a spherically symmetric trap with incomensurable frequencies $`\omega ^g`$ and $`\omega ^e`$, for the ground and the excited state respectively. The latter assumption simplifies the dynamics of the spontaneous emission processes in the FL limit. The cooling process consists of sequences of Raman pulses of frequencies, adjusted in such a way that they induce the transition of atoms to the lower motional states of the trap. We assume that the power of cooling lasers is sufficiently weak, and therefore during each pulse no significant population in $`|e`$ is present. This allows to eliminate adiabatically the level $`|e`$, and to consider only the density matrix $`\rho (t)`$ for the atoms in $`|g`$, and being diagonal in the Fock representation corresponding to the bare trap levels.

We assume that the laser acts only on one component, whereas the other one is cooled sympathetically, which is sufficient to reach $`TT_F`$. Using the standard theory of quantum-stochastic processes GardinerB ; Carmichael  we derive the quantum master equation (ME) for the density matrix $`\rho (t)`$ in a similar way as that for bosons CoolBosons :

$$\dot{\rho }(t)=_0\rho +_1\rho +_2\rho ,$$
(1)
where $`_0\rho =i\widehat{H}_{\mathrm{eff}}\rho (t)+i\rho (t)\widehat{H}_{\mathrm{eff}}^{}+𝒥\rho (t)`$, $`_1\rho =i[\widehat{H}_{\mathrm{las}},\rho (t)]`$, $`_2\rho =i[\widehat{H}_{\mathrm{coll}},\rho (t)]`$. The first term in the ME (1) describes the evolution of the atoms in the vacuum of the electromagnetic field. The effective non-Hermitian Hamiltonian $`H_{\mathrm{eff}}`$ is defined as

$`\widehat{H}_{\mathrm{eff}}`$ $`=`$ $`{\displaystyle \underset{𝐦}{}}\omega _𝐦^gg_𝐦^{}g_𝐦+{\displaystyle \underset{𝐥}{}}(\omega _𝐥^e\delta )e_𝐥^{}e_𝐥+`$ (2)
$`+{\displaystyle \underset{𝐧}{}}\omega _𝐧^bb_𝐧^{}b_𝐧i{\displaystyle \frac{\gamma }{2}}{\displaystyle \underset{𝐥,𝐦}{}}\xi _{\mathrm{𝐥𝐦}}e_𝐥^{}g_𝐦g_𝐦^{}e_𝐥,`$
and the so-called jump super-operator $`𝒥`$ is given by

$$𝒥\rho (t)=\underset{𝐥,𝐦}{}\xi _{\mathrm{𝐥𝐦}}g_𝐦^{}e_𝐥\rho (t)e_𝐥^{}g_𝐦.$$
(3)
Here $`g_𝐦`$ ($`e_𝐥`$) is the annihilation operator of atoms of the first component, in the ground (excited) state, and in the trap level $`𝐦`$ ($`𝐥`$). The annihilation operator of atoms of the second component, and in the trap level $`𝐧`$, is denoted by $`b_𝐧`$. These operators fulfill the standard fermionic anticommutation relations: $`\{g_𝐦,g_𝐧^{}\}=\{e_𝐦,e_𝐧^{}\}=\{b_𝐦,b_𝐧^{}\}=\delta _{𝐦,𝐧}`$. In Eq. (2) $`\gamma `$ denotes the single-atom effective spontaneous emission rate, $`\omega _𝐦^g`$, $`\omega _𝐦^e`$, $`\omega _𝐦^b`$ are the energies of state $`𝐦`$ of the trap: $`|g`$, $`|e`$, $`|b`$ respectively, and $`\delta `$ is the laser detuning from the atomic transition. The coefficients $`\xi _{\mathrm{𝐥𝐦}}`$ are defined as follows $`\xi _{\mathrm{𝐥𝐦}}=_0^{2\pi }𝑑\varphi _0^\pi 𝑑\theta \mathrm{sin}\theta 𝒲(\theta ,\varphi )|\eta _{\mathrm{𝐥𝐦}}(\stackrel{}{k})|^2`$, where $`𝒲(\theta ,\varphi )`$ is the fluorescence dipole pattern, and $`\eta _{lm}(\stackrel{}{k})=e,𝐥|e^{i\stackrel{}{k}\stackrel{}{r}}|g,𝐦`$ the Frank-Condon factors.

The interaction of the laser light with the atoms is governed by the Hamiltonian $`\widehat{H}_{\mathrm{las}}`$:

$$\widehat{H}_{\mathrm{las}}=\frac{\mathrm{\Omega }}{2}\underset{𝐥,𝐦}{}\eta _{\mathrm{𝐥𝐦}}(\stackrel{}{k}_L)e_𝐥^{}g_𝐦+H.c.,$$
(4)
where $`\mathrm{\Omega }`$ is the Rabi frequency and $`\stackrel{}{k}_L`$ is the wavevector of the cooling laser. Binary collisions are described by

$$\widehat{H}_{\mathrm{coll}}=\underset{𝐦,𝐧,𝐪,𝐩}{}U_{𝐦,𝐧,𝐪,𝐩}g_𝐦^{}g_𝐧b_𝐪^{}b_𝐩,$$
(5)
where due to Fermi statistics only collisions between different species are allowed. We neglect the collisions with the atoms in $`|e`$, since only a small fraction of atoms is excited in each pulse. The collisional amplitudes are defined as $`U_{𝐦,𝐧,𝐪,𝐩}=(4\pi \mathrm{}a_{\mathrm{sc}}/m)_{R^3}d^3x\varphi _𝐦^{}(𝐱)\varphi _𝐧(𝐱)\beta _𝐪^{}(𝐱)\beta _𝐩(𝐱)`$, where $`\varphi _𝐧`$ ($`\beta _𝐧`$) denotes the wavefunction of the state $`𝐧`$ in the $`|g`$, ($`|b`$) trap, and $`a_{\mathrm{sc}}`$ is the scattering length.

We require the FL condition $`\gamma <\omega `$ at least in the initial phase of the cooling process, and also assume that $`\mathrm{\Omega }<\gamma `$ and $`\mathrm{\Omega }^2/\gamma \omega `$, which allows for the adiabatic elimination of $`|e`$, performed by means of projection operator techniques GardinerB . For weak interactions between the two species, the dynamics due to the collisions and the dynamics due to the laser-cooling are independent Coll . The laser cooling may be described by the ME without the collisional part, whereas the collisional part takes the form of a quantum Boltzmann master equation (QBME) QBME . The probability of a laser-induced transition from state $`𝐦`$ to state $`𝐧`$ of the trap $`|g`$ is given by

$`P_{𝐦𝐧}^{\mathrm{opt}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }^2}{\gamma }}N_𝐦^1(1N_𝐧^1)\times `$ (6)
$`\times {\displaystyle \underset{l}{}}{\displaystyle \frac{\gamma ^2\xi _{\mathrm{𝐥𝐧}}|\eta _{\mathrm{𝐥𝐦}}(\stackrel{}{k}_L)|^2}{4[\delta (\omega _𝐥^e\omega _𝐦^g)]^2+\gamma ^2R_{\mathrm{𝐦𝐥}}^{}{}_{}{}^{2}}},`$
where $`R_{\mathrm{𝐦𝐥}}=_𝐧\xi _{\mathrm{𝐥𝐧}}(1N_𝐧^1+\delta _{𝐧,𝐦})`$ are the factors modifying the spontaneous emission rate in a Fermi gas, and $`N_𝐦^1`$ and $`N_𝐦^2`$ denotes the number of atoms occupying the state $`𝐦`$ of the trap $`|g`$ and $`|b`$, respectively. In the regime of quantum degeneracy, coefficients $`R_{\mathrm{𝐦𝐥}}`$ vanish, inhibiting the spontaneous emission and forcing the atoms to remain excited for a long time. This prolongs the cooling process, and can result in excited-ground collisions, which lead to heating and losses. In addition, the adiabatic elimination used in the derivation of Eq. (6) ceases to be valid. The cooling efficiency is also decreased due to the fermionic inhibition factor $`(1N_𝐧)`$ in the numerator of (6). The negative influence of the statistics can be overcome by dynamically increasing $`\gamma `$ during the Raman cooling fermicool , in order to avoid the inhibition effects, but still remaining in the FL regime: $`\gamma R_{\mathrm{𝐦𝐥}}<\omega `$. Still, a small fraction of the atoms will remain in the excited state after the cooling pulses, and has to be removed from the trap in order to avoid inelastic collisions. The latter aim can be achieved by optically pumping the excited atoms to a third non-trapped level.

The probability of a collision between two fermions of different species, from the states $`𝐧`$ and $`𝐩`$ to the states $`𝐦`$ and $`𝐪`$, respectively, is given by

$`P_{𝐧,𝐩𝐦,𝐪}^{\mathrm{coll}}`$ $`=`$ $`{\displaystyle \frac{\pi }{\omega }}N_𝐧^1N_𝐩^2(1N_𝐦^1)(1N_𝐪^2)\times `$ (7)
$`\times \left|U_{𝐦,𝐧,𝐪,𝐩}\right|^2\delta _{E_𝐧+E_𝐩,E_𝐦+E_𝐪}.`$
The fermionic inhibition factors $`(1N_𝐦^1)`$ also slow the collisional processes.

In our simulations we assumed the atoms as confined in an isotropic, 3D harmonic trap with frequency $`\omega =\omega ^g=\omega ^e=\omega ^b`$. Due to the limitations when simulating relatively large systems ($`N10^4`$), we employ ergodic approximation, i.e. we assume that the populations of the states with the same energy are equal. This approximation relies on the fact that thermalization inside the same energy shell is much faster than between different energy shells. In addition we assume that the collisional processes is much faster than the laser cooling. The detailed calculation of the transition probabilities within ergodic approximation is presented in Ref. fermicool2 . Due to numerical limitations we assumed $`\eta =2`$, which could be e.g. the case of potassium atoms in a dipole trap with $`\omega =2\pi \times 2.4`$ kHz (employed in Li experiments Granade ), and a laser wavelength $`\lambda 720`$ nm. We consider a scattering length for the interactions between the two species $`a_{\mathrm{sc}}=157a_0`$, where $`a_0`$ is the Bohr radius. This value corresponds to the interactions between $`|F=9/2,m_F=9/2`$ and $`|F=9/2,m_F=7/2`$ of 40K DeMarco . The trap is assumed to have the same depth for both species, and contains $`81`$ energy levels ($`91881`$ states). We assume both components to have equal initial number of atoms, $`N=10660`$, corresponding to a Fermi energy $`E_F=38\mathrm{}\omega `$. For this relatively large number of atoms, the dynamics generated by the collisional part of the ME equation leads to an equilibrium distribution, which agrees very well with the one calculated from the grand canonical ensemble fermicool2 . Hence we start the simulations from a thermal distribution.

In a Fermi gas the main loss sources are provided by background collisions, resulting from non ideal vacuum conditions, and photoassociation. The photoassociation losses (when the laser is tuned between molecular resonances) are typically of the order of $`10^{14}`$cm3/s for laser intensities of $`1`$mW/cm2 Machholm . In our case the applied laser intensities are typically $`1000`$ times smaller and we estimate that for $`N=10660`$ atoms, the atomic density is smaller than $`3.5\times 10^{14}`$cm-3. Therefore the photoassociation losses can be safely neglected. On the contrary, the background collision apart from decreasing the number of atoms, generate holes deep within the Fermi sea. Those holes, after subsequent thermalization, may lead to a significant heating Timmer . We have assumed that background losses depopulate each state of the trap with the rate $`\gamma _{\mathrm{bg}}`$, which is independent of the energy of the state

$$\dot{N}_n^j=\gamma _{\mathrm{bg}}N_n^j,$$
(8)
where $`j=1,2`$ enumerates the components. In the simulations we assume $`\gamma _{\mathrm{bg}}=1/350`$ Hz HollandF . We include also in the simulations the losses produced by the removal of long-living excited atoms at the end of each cooling pulse.

Fig. 1 shows the evolution of the temperature of the laser-cooled two-component Fermi gas. Initially we consider $`T_0=T_F`$. The cooling process was divided into four stages, each one consisting of a sequence of two Raman pulses. The employed pulses are characterized by the following parameters: detuning $`\delta /\omega =\{(11,12)`$, $`(16,17)`$, $`(19,20)`$,$`(25,26)\}`$ respectively, Rabi frequency $`\mathrm{\Omega }/\gamma =\{(0.113,0.113)`$, $`(0.008,0.012)`$, $`(0.0025,0.004)`$,$`(0.003,0.01)\}`$ respectively, and length $`\mathrm{\Delta }t/\omega ^1=\{(250,250)`$, $`(2000,2000)`$, $`(4000,4000)`$,$`(4000,6500)\}`$ respectively. For these parameters not more than $`10\%`$ of the atoms is excited during each pulse, and thus the conditions of the adiabatic elimination are fulfilled. The temperature was determined by fitting the calculated distribution of fermions to a thermal distribution. As one can observe, a final temperature $`T0.008T_F`$ may be reached within $`8`$ s. The losses associated with the cooling process do not exceed $`2`$% and are slightly larger for the laser-cooled component, due to the removal of the long-living excited atoms (inset of Fig. 1).

We have analyzed the effect of the losses on the laser-cooled gas. To this aim we have considered an initial gas at $`T=0.008T_F`$ (end of the cooling process in Fig. 1), and compared two different cases: (i) the laser is turned off, and the gas is heated due to the background collisions, (ii) the laser is turned on, and the cooling pulses of the last stage are continuously applied. Fig. 2 presents the evolution of the temperature of both species for these two cases. As observed, the laser cooling compensates for the heating induced by the creation of holes in the degenerate distribution. Hence, it helps to maintain the degenerate gas for a relatively long time in the trap. We have verified that the continuous application of the laser cooling does not lead to substantially larger losses.

Finally, we have analyzed the results for the heating due to background collisions with the help of an analytical model. For the background losses that decrease the population of the trap levels in the way described by Eq. (8), one can calculate the evolution of the temperature analytically Timmer . To this end, it is sufficient to observe that the Eq. (8) do not change the mean energy per particle $`\epsilon `$, which has to remain constant in the absence of cooling: $`\epsilon (N(t),T(t))=`$ const. At low temperatures $`\stackrel{~}{T}=T/T_F1`$, the mean energy per particle in a harmonic trap is given by $`\epsilon (N,T)=(3/4)E_F(N)+(\pi ^2/2)(k_BT)^2/E_F(N)`$. Using $`\epsilon (N(t),T(t))=\epsilon (N_0,T_0)`$, where $`N_0`$ is the initial number of atoms, and $`T_0`$ is the initial temperature of the gas, we obtain

$$\stackrel{~}{T}(t)=\sqrt{\frac{E_F(N_0)}{E_F(N(t))}\stackrel{~}{T}_0^2+\frac{3}{2\pi ^2}\left(\frac{E_F(N_0)}{E_F(N(t))}1\right)}.$$
(9)
For a large system $`E_F(N)\mathrm{}\omega (6N)^{1/3}`$, and $`N(t)N_0\mathrm{exp}(\gamma _{\mathrm{bg}}t)`$, hence

$$\stackrel{~}{T}(t)=\sqrt{e^{\gamma _{\mathrm{bg}}t/3}\stackrel{~}{T}_0^2+3\left(e^{\gamma _{\mathrm{bg}}t/3}1\right)/2\pi ^2}.$$
(10)
Fig. 3 presents the evolution of the temperature for a rate $`\gamma _{\mathrm{bg}}=10`$ Hz, which is much larger than the typical experimental rates, and was employed to reduce the amount of time-consuming numerical calculations. Fig. 3 compares our numerical results, the predictions of the analytical formula of Timmer  derived for a homogeneous system, and the data calculated from Eq. (9) with $`E_F(N)\mathrm{}\omega \left((6N)^{1/3}2\right)`$, which is more appropriate for a finite size system TimmHarm . The analytic curve fits very well to the numerical data. The small discrepancy for larger temperatures originate from the fact, that the expression for the mean energy, which we use, is valid up to order $`(T/T_F)^2`$. From the plot, it is clearly seen that the increase of the temperature in a trapped gas is much slower that in the homogeneous case.
In conclusion, we have studied the laser cooling of trapped two-component Fermi gases. Our method exploits the dynamical adjusting of the effective spontaneous emission rate in Raman cooling. Using Monte Carlo simulations, we have shown that laser cooling is able to cool the two-component Fermi system below $`0.01`$ $`T_F`$. We have also discussed the losses which may affect the laser-cooled gas. In this context, we have shown that our laser-cooling scheme can be employed to maintain a two-component Fermi gas at a fixed temperature in the presence of background collisions, without significant additional atom losses. Finally we have derived an analytic formula for the temperature of a trapped Fermi gas heated by background collisions. We have shown that the heating rate is smaller than in the homogeneous case.

We acknowledge support from the Alexander von Humboldt Stiftung, the Deutsche Forschungsgemeinschaft, the RTN Cold Quantum gases, the ESF Program BEC2000+, the Polish KBN Grant No. 5-P03B-103-20, and the Russian Foundation for Fundamental Research.


text:
# Diffusion and viscosity in a supercooled polydisperse system

## I Introduction

The rapid cooling of a liquid below its freezing temperature transforms it into a long-lived metastable amorphous solid or glass. Understanding the dynamics of the system near the glass transition is an intense field of research since the last few decades. There have been many experimental as well as simulation studies which focus on the dynamics of dense supercooled liquids well above the glass transition and also near the glass transition temperature ($`T_g`$). The basic aim of all these studies was to characterise quantitatively the observed very complex dynamics of the system as it approaches the glass transition from the above. Close to the glass transition, the shear viscosity ($`\eta `$) and the microscopic structural relaxation time ($`\tau `$) of the so-called fragile glass-forming liquids show divergence with a strongly non-Arrhenius temperature dependence. This divergence is often well represented by the Vogel-Fulcher-Tammann (VFT) equation

$$\eta (T)=A_\eta \mathrm{exp}[E_\eta /(TT_o^\eta )]$$
(1)
where $`A_\eta `$ and $`E_\eta `$ are temperature-independent constants and $`T_o^\eta `$ ($`<T_g`$) is the temperature at which $`\eta `$ diverges. Note that at low temperatures the increasingly slow dynamics of the so-called fragile liquids is simultaneously manifested in the stretched exponential decay of the stress correlation function (with a strongly temperature dependent stretching parameter). The VFT dependence (Eq. 1) is thus accompanied by the strong non-exponential relaxation observed near the glass transition temperature.

The dramatic slow down of the dynamics near the glass transition is not well understood and still remains the most challenging problem in the physics of glasses. Several theories have been proposed to understand the anomalous relaxation dynamics of deeply supercooled liquids. Although the ideal version of the non-linear mode coupling theory (MCT) gives a microscopic picture of this slowing down, it predicts a structural arrest, i.e., a transition from ergodic to nonergodic behavior, at a critical temperature $`T_c`$, well above the laboratory glass transition temperature $`T_g`$. Near $`T_c`$, the importance of the influence of potential energy landscape on the relaxation processes has now widely been accepted and the strongly correlated jump motion is observed to be the dominant mode for mass transport, which is not included in ideal MCT.

Another important characteristic feature in the dynamics of deeply supercooled liquids is the decoupling between translation diffusion and the shear viscosity of the medium. At high temperature, over a wide range of liquid states, the translational diffusion is inversely proportional to viscosity, in accordance to the Stokes-Einstein (SE) relation given by

$$D_T=\frac{k_BT}{C\pi \eta R}$$
(2)
where $`R`$ is the spherical radius of the diffusing particle and $`C`$ is a numerical constant that depends on the hydrodynamic boundary condition. However, several recent experimental and simulation studies on strongly supercooled ’fragile’ glass forming liquids have shown significant deviations from the SE relation. As the temperature is lowered toward $`T_g`$, it is found that the translational diffusion is larger than the value predicted by the SE relation and in some cases even two to three orders of magnitude larger. The enhanced diffusion at low temperatures is sometimes explained in terms of a power law behavior $`D_T\eta ^\alpha `$ with $`\alpha <1`$. Both the experiments and simulation studies have recently been evidenced the enhancement of the translational diffusion coefficient is due to the spatially heterogeneous dynamics in deeply supercooled liquids.

Computer simulations have played a key role in augmenting our understanding of various aspects of the dynamics of the supercooled liquids from a microscopic viewpoint. Unfortunately, the simple one-component systems such as soft or hard spheres or Lennard-Jones systems crystallize rapidly on lowering the temperature below the melting point ($`T_m`$) and, therefore, cannot be utilized as a model for studying the complex dynamical behavior near the glass transition temperature. A natural way to avoid crystallization is to use binary mixtures of atoms with different diameters. A large number of molecular dynamics (MD) simulations have recently been carried out in supercooled model binary mixtures near the glass transition as well as below the glass transition temperature.

However, one is often interested in the consequences of the disorder introduced by the dissimilarity of the particles. Synthetic colloids, by their very nature, frequently exhibit considerable size polydispersity. Polydispersity is also common in industrially produced polymers which always contain macromolecules with a rang of chain length. Colloidal particles are an excellent model of hard spheres and perhaps the simplest possible experimental system of interacting particles to study the glass transition. Several experiments and simulations have shown that the crystal phase of the colloidal systems can exist as a thermodynamically stable phase only for polydispersities (standard deviation of the size distribution divided by the mean) less than a ’terminal’ value, in the range of 0.05-0.15.

Interestingly, recent experiments on colloidal supercooled fluids and colloidal glasses allowed one to obtain the information on the microscopic details of the dynamics of the individual particles. These experiments have shown the presence of dynamic heterogeneity in deeply supercooled colloidal systems. The motion of the relatively fast moving particles is found to be highly correlated and form connected clusters, whose size increases as one approach the glass transition. More recently, Sear has carried out MD simulation of a dense polydisperse hard sphere fluid to study the effect of polydispersity on the slow dynamics. The simulation results also show the clustering of the fast-moving particles in agreement with the experiments, though the dynamics appears to be less heterogeneous. The heterogeneous nature of the dynamics has also been observed in the Monte Carlo simulation study of polydisperse hard spheres close to the glass transition.

It is worth mentioning that besides the short-range hard-core interaction, addition of a nonadsorbing soluble polymer in the stable colloidal suspension gives rise to a weak, long-range attraction between the colloidal particles by means of the depletion interaction. Recently, the consequences of this attractive interaction on the glass transition are nicely explained in a combined experimental, theoretical, and simulation study by Pham et.al. Interestingly, with increase in the strength of the short-range attractive interaction, two qualitatively different glassy states are found with a reentrant glass transition line.

The size distribution in real colloids generally leads to a distribution in mass of the particles. The importance of the mass polydispersity on the dynamics of a realistic system having size polydispersity has recently been analyzed in molecular dynamics simulation study of a Lennard-Jones (LJ) polydisperse fluid near the triple point of the corresponding monodisperse LJ system. Polydispersity is commonly found in many systems of industrial applications and to mimic the interparticle interactions, the Lennard-Jones potential generally serve as a good starting model. Thus, it will be a general interest to study the impact of polydispersity on the dynamics of a deeply supercooled polydisperse fluid, where particles interacting via the LJ potential. More importantly, this will enable us to compare the properties of the system with model binary LJ mixtures whose dynamics near the glass transition has been studied extensively in simulations.

In this study, we have performed extensive molecular dynamics simulations of a system of polydisperse LJ spheres with continuous range of diameters and mass. The temperature dependence of the dynamic properties such as viscosity ($`\eta `$) and the self-diffusion coefficients ($`D_i`$) for the different size particles is studied by varying the temperature ($`T`$) over a large range at a constant high pressure ($`P`$). Both the viscosity and diffusion show super-Arrhenius temperature dependence and the calculated value of the fragility parameter (D) show that the present system is more fragile than the well-known Kob-Andersen binary mixture. The critical temperature obtained from the VFT fit to the diffusions ($`T_o^{D_i}`$) show strong dependency on the radius ($`R`$) of the particles. In addition, the critical temperature obtained from the VFT fit to the viscosity ($`T_o^\eta `$) is much higher when compared to those of diffusion coefficients, where the deviation is largest for the smaller size particles. This clearly reflects the deviation from the Stokesian diffusion in the proximity of the glass transition temperature. Most interestingly, at lower temperatures the diffusion shows a highly nonlinear size dependence when plotted against the inverse of the radius ($`R`$) of the particles. The reason for the breakdown of the Stokes-Einstein (SE) relation can be analyzed from the trajectories of the particles. We find that at low temperature the hopping processes being the primary mode of particle diffusion for both the smaller and bigger size particles.

The organization of the rest of the paper is as follows. In the next section, we describe in detail the system studied here and the details of the simulations. The simulation results are analyzed and discussed in section III. Finally, we end with some concluding remarks in section IV.

## II System and Simulation Details

We have performed a series of equilibrium isothermal-isobaric ensemble (N P T) molecular dynamics simulations in three dimensions of a system of $`N=256`$ particles of mean radius $`\overline{\sigma }`$ and mass $`\overline{m}`$ with polydispersity in both size and mass. The interaction between any two particles is modeled by means of shifted force Lennard-Jones (LJ) pair potential, where the standard LJ is given by

$$u_{ij}^{LJ}=4ϵ_{ij}\left[\left(\frac{\sigma _{ij}}{r_{ij}}\right)^{12}\left(\frac{\sigma _{ij}}{r_{ij}}\right)^6\right]$$
(3)
where $`i`$ and $`j`$ denotes two different particles and $`\sigma _{ij}=(\sigma _i+\sigma _j)/2`$ where $`\sigma _i`$, $`\sigma _j`$ are the diameters of the particles $`i`$ and $`j`$, respectively. In the shifted-force potential both the potential and force are continuous at a cutoff radius $`r_c`$ and we choose a value of $`r_c=2.5\overline{\sigma }`$.

The polydispersity in size is introduced by random sampling from the Gaussian distribution of particle diameters $`\sigma `$

$$P(\sigma )=\frac{1}{\delta \sqrt{2\pi }}\mathrm{exp}\left[\frac{1}{2}\left(\frac{\sigma \overline{\sigma }}{\delta }\right)^2\right],$$
(4)
where $`\delta `$ is the width of the distribution. The standard deviation ($`\delta `$) of the distribution divided by its mean $`\overline{\sigma }`$ gives a dimensionless parameter, the polydispersity index $`s=\delta /\overline{\sigma }`$. The simulations are carried out here for a fixed value of the polydispersity index $`s=0.1`$. The masses of the particles are varied accordingly and we assume that the mass of a particle $`i`$ is scaled by its diameter as $`m_i=\overline{m}(\sigma _i/\overline{\sigma })^3`$. The LJ energy parameter $`ϵ_{ij}`$ is assumed to be same for all particle pairs and denoted as $`ϵ`$. All the quantities in this study are given in reduced units, that is, length in units of $`\overline{\sigma }`$, temperature $`T`$ in units of $`ϵ/k_B`$, pressure $`P`$ in units of $`ϵ/\overline{\sigma }^3`$, and the unit of time is $`\tau =\sqrt{\overline{m}\overline{\sigma }^2/ϵ}`$. Note that if one assumes the argon units then $`\tau `$ = 2.2 ps.

All simulations in the NPT ensemble were performed using the Nose-Hoover-Andersen method, where the external reduced temperature ($`T^{}`$) is varied over a large range from 1.3 to 0.67 keeping the external reduced pressure ($`P^{}`$) fixed at 10.0. Throughout the course of the simulations, the barostat and system’s degrees of freedom are coupled to an independent Nose-Hoover chain (NHC) of thermostats, each of length 5. The extended system equations of motion are integrated using the reversible integrator method. The higher order multiple time step method has been employed in the NHC evolution operator which lead to stable energy conservation for non-Hamiltonian dynamical systems. The extended system time scale parameter used in the calculations was taken to be $`0.93`$ for $`T^{}1.0`$ and $`1.16`$ for $`T^{}<1.0`$ for both the barostat and thermostats.

A time step of $`0.001`$ is employed for $`T^{}1.0`$ and $`0.002`$ for $`T^{}<1.0`$. The equilibration and data collection steps are also varied accordingly depending upon the temperature of the system. For $`T^{}1.0`$, the equilibration steps are varied from $`5\times 10^5`$ to $`10^6`$ and the data collection steps are $`10^6`$, whereas for $`T^{}<1.0`$, the equilibration steps are varied from $`5.0\times 10^5`$ to $`2.0\times 10^6`$ and the data collection steps from $`10^6`$ to $`2.5\times 10^7`$. At each temperature, all the dynamic quantities are averaged over five independent runs. Diffusion coefficients ($`D_i`$) for the different size particles are calculated from the slope of the corresponding mean square displacements (MSD) in the diffusive limit and viscosity is calculated from the auto-correlation of the off-diagonal components of the microscopic stress tensor, via the standard Green-Kubo formula. As the system is isotropic, we have taken an average over three different off-diagonal stress correlations for each of the five data set.

## III Results and Discussion

In order to make sure that there is no crystallization, we have calculated the radial distribution functions $`g(r)`$ which describes the average structure of the fluid. The radial distribution function calculated for $`T^{}=0.67`$ is shown in figure 1. The decay of correlations with increase in distance clearly reflects the absence of any long range order, a characteristic feature of the fluid.

The plot of ln($`\eta `$) as a function of inverse of temperature ($`1/T^{}`$) in figure 2(a) clearly shows a super-Arrhenius behavior of the viscosity. In figure 2(b) we show a VFT fit to the viscosity (Eq. 1) by plotting ln($`\eta `$) against (1/($`TT_o^\eta `$)) where $`T_o^\eta `$ is equal to 0.57. As in other fragile liquids, it shows that the divergence of viscosity is quite well described by the VFT equation. From the fitting we obtain the values of $`A_\eta `$ and $`E_\eta `$ as 2.0 and 0.81, respectively. We have also calculated the fragility parameter (D = $`E_\eta /T_o^\eta `$) as defined by Angell. Using the values of the fitting parameters ($`E_\eta `$ and $`T_o^\eta `$) obtained within the temperature range investigated, its value is $`1.42`$. This classifies the present system into a strongly fragile liquid. Thus, the dense random packing of unequal size particles makes the present system more fragile when compared with a recent simulation study on Kob-Andersen binary mixture (D $`2.45`$).

In a polydisperse system, all the particles are unequal in size, so their diffusion coefficients will also differ. We categorize the particles into different subsets where particles of diameters within $`0.05\overline{\sigma }`$ are assumed to be members of the same subset. For the polydispersity index $`s=0.1`$, we find that the minimum and maximum diameter of the particles are $`0.75\overline{\sigma }`$ and $`1.25\overline{\sigma }`$, respectively. Thus subsets of particles with diameters in the ranges $`0.75`$ to $`0.8\overline{\sigma }`$ and $`1.2`$ to $`1.25\overline{\sigma }`$ corresponds to smallest and largest spheres, respectively. The diffusion coefficients for different subsets of particles are calculated at each temperature. It is well-known that in deeply supercooled liquids the diffusion coefficient shows non-Arrhenius temperature dependence and can be fitted by a VFT law

$$D_i(T)=A_{D_i}\mathrm{exp}[E_{D_i}/(TT_o^{D_i})],$$
(5)
where the index $`i`$ stands for the different subsets of particles. $`T_o^{D_i}`$ is the critical temperatures for $`i`$-th species at which the diffusion coefficients ($`D_i`$) vanishes. The diffusion coefficients for each subsets of particles have been fitted to the above equation and we show the plot of ln($`D_i`$) against $`1/(TT_o^{D_i})`$, for the smallest ($`D_1`$) and largest spheres ($`D_{10}`$) in figures 3(a) and 3(b), respectively, where $`T_o^{D_1}`$ and $`T_o^{D_{10}}`$ are 0.46 and 0.5. Thus the critical temperatures ($`T_o^{D_i}`$) depends on the size of the particles and increases with size of the particles $`\overline{R}_i`$ ($`\overline{R}_i`$ is the mean radius of the $`i`$-th subset). This is shown in figure 4. For the largest spheres, the critical temperature ($`T_o^{D_{10}}`$) is smaller than the corresponding critical temperature obtained from VFT fit to the viscosity ($`T_o^\eta =0.57`$). This clearly signifies that near the glass transition the diffusion is partly decoupled from the viscosity and for smaller particles the degree of decoupling is more. The smaller particles remain mobile even when bigger particles are almost frozen.

In figure 5 we plot the diffusion constants ($`D_i`$) against $`1/\overline{R}_i`$ at the lowest temperature of $`T^{}=0.67`$ and compared with the well-known hydrodynamic Stokes-Einstein (SE) relation (Eq. 2 with $`C`$ = 6, the stick boundary condition). It clearly shows the markedly non-Stokesian behavior of the diffusion at low temperatures. Interestingly, the fitting to the simulated data points show a highly nonlinear size dependence of the diffusion. This is a clear evidence that the breakdown of SE law is more severe for the smaller size particles. In order to get an estimate of the degree of decoupling (between diffusion and viscosity) for the smallest size particles, we have fit the inverse diffusion coefficient ($`1/D_1`$) versus $`\eta /T`$. While at high T, it asymptotically satisfies the SE relation (slope is 1), the fit to the low temperature data gives the slope $`\alpha 0.5`$ (that is, diffusion shows the power law behavior $`D_1\eta ^{0.5}`$). It should be noted that the dynamics of a polydisperse liquid is more heterogeneous than a monodisperse or bidisperse system due to the different time scales involved for different sizes and masses of the particles. The smaller particles are on average faster than others over all time scales. This becomes more prominent at lower temperatures where the relaxation time of the system is very high. At low temperature, the observed nonlinear dependence of diffusion on size is related to the increase in dynamic heterogeneity in a polydisperse system.

A more detailed analysis of the diffusion can be obtained from a closer examination of the self-part of the van Hove correlation function $`G_s(r,t)`$. This gives the distribution of the displacements ($`r`$) of a particle in a time interval $`t`$. We calculate $`G_s(r,t)`$ for the smallest ($`\sigma `$ = 0.75 to 0.8$`\overline{\sigma }`$) and largest ($`\sigma `$ = 1.2 to 1.25$`\overline{\sigma }`$) particles for different time intervals at $`T^{}=0.67`$, the lowest temperature investigated. Figures 6(a) and 6(b) shows the correlations for the smallest and largest spheres, respectively. For the smallest particles, a gradual development of a well-defined second peak at $`r1.0\overline{\sigma }`$ is clearly visible with increase in time (figure 6(a)). However, for the largest particles, the distribution becomes bimodal at relatively longer time scales (figure 6(b)). The occurrence of the secondary peak, observed also in other model binary mixtures at low temperatures, is an evidence of the jump motion in the dynamics of the particles.

To characterize the single particle dynamics further, we have evaluated the self-intermediate scattering function $`F_s(k,t)`$, the spatial Fourier transform of $`G_s(r,t)`$, for different subsets of the particles for a fixed value of the reduced wave number $`k\overline{\sigma }2\pi `$ at $`T^{}=0.67`$. The long time decay of $`F_s(k,t)`$ is well fitted by the Kohlrausch-Williams-Watts (KWW) stretched exponential form

$$F_s^i(k,t)=\mathrm{exp}(\frac{t}{\tau _i})^{\beta _i}$$
(6)
where $`\tau _i`$ and $`\beta _i`$ are the relaxation time and the stretching exponent of the $`i`$-th subset. We find that both $`\tau _i`$ and $`\beta _i`$ increases with an increase in the size of the particles, as has been observed earlier by other authors in binary mixtures. The $`F_s(k,t)`$ calculated for the smallest (subset 1) and largest (subset 10) particles along with the KWW fits is shown in figures 7(a) and 7(b), respectively. Note that we fit the functions $`[F_s^i(k,tt_o)/F_s^i(k,t_o)]`$ ($`t>t_o`$), to the KWW form to quantify their long time behavior. For the smallest particles, the values of the fitting parameters are found to be $`\tau _1242`$ and $`\beta _10.49`$, whereas for the largest particles they are $`\tau _{10}717`$ and $`\beta _{10}0.64`$. The enhanced stretching ($`\beta _1<\beta _{10}`$) at long times is due to the more heterogeneity probed by the smaller size particles than that by the larger size particles during the time scale of decay of their $`F_s(k,t)`$.

In order to determine the extent of the jump-type motion more clearly we follow the trajectory of the individual particles. A close inspection of the simulated trajectory of the smallest and largest particles reveals several interesting features. Figure 8(a) display the projections onto an x-y plane of the trajectory of a typical smallest size particle over a time interval $`\mathrm{\Delta }t=500\tau `$ and figure 8(b) shows the paths followed by a largest size particle over time interval $`\mathrm{\Delta }t=2000\tau `$, both at $`T^{}=0.67`$. At this temperature the dynamics is dominated by ’hopping’; particles remain trapped in transient cages created by the surrounding particles for quite some time and then moves significant distance (approximately one interparticle distance) by making a jump to a new cage. For the larger size particles, the jump motion begins to take place at relatively longer time scales (compare also figures 6(a) and 6(b)). Thus, it clearly elucidates that in a system with particles of all different sizes and masses, the hopping is the dominant diffusive mode both for the smaller and bigger size particles. It is the frequent hopping in case of smaller size particles that leads to the severe breakdown of SE relation.

## IV Conclusions

In summary, we have presented the results of large scale computer simulations for a supercooled polydisperse system with large variations in temperature at a fixed high pressure. The interparticle interaction is represented by the standard Lennard-Jones (LJ) potential and the particles in the system are considered to be polydisperse both in size and mass. Characteristic of a fragile glass former, the super-Arrhenius temperature dependence is observed for the viscosity and also for the self-diffusion coefficients of different size particles. Furthermore, within the temperature range investigated, the value obtained for the Angell’s fragility parameter (D $`1.4`$) establish the present system as a strongly fragile liquid.

In a dense polydisperse system, a wide range of time scales are involved due to the different size and mass of the particles. Thus, upon lowering the temperature, the dynamics is increasingly more heterogeneous compared to a monodisperse or bidisperse system. The increase of critical temperature for diffusion ($`T_o^{D_i}`$) with the size of the particles suggests that the dynamics is indeed heterogeneous. In addition, the critical temperature for viscosity ($`T_o^\eta `$) is found to be larger than that for the diffusion, indicating the decoupling of translational diffusion from the viscosity commonly observed in deeply supercooled liquids. Interestingly, a marked deviation from the Stokesian diffusion is observed where the dependence on size of the particles is highly nonlinear.

The analysis of the self-part of the van Hove correlation functions $`G_s(r,t)`$ showed a clear signature of the single particle ’hopping’ (indicated by the second neighbor peak) in the dynamics of both the smallest and largest size particles at low temperatures. However, for the larger particles, the hopping processes are found to occur at relatively longer times. The relevance of the hopping processes at low temperatures is further investigated in detail by following the trajectory of the individual particles. Crossover from continuous Brownian to hopping motion takes place at shorter time scales for the smaller size particles than that for the larger size particles.

In the present system the size of all the particles are different. It would be interesting to see whether the jump motion executed by the individual particles occurs over a single energy barrier or it takes place via a number of ”intermediate” inherent structures in the potential energy landscape. A recent molecular dynamics simulations on LJ binary mixture have shown that such a transition does not correspond to transitions of the system over single energy barriers. In addition, in a deeply supercooled liquid the jump motions are associated with strong nearest-neighbour correlations, in which several neighboring atoms jump at successive close times. It is to be noted that similar correlations have been observed here also. Recently, a computer simulation study of a deeply supercooled binary mixture has shown that the local anisotropy in the stress is responsible (at least partly) for the particle hopping. However, the molecular origin of the jump motions observed here (highly disordered system) is not clear and we are presently pursuing this problem.

Acknowledgments

This work was supported in part by the Council of Scientific and Industrial Research (CSIR), India and the Department of atomic energy (DAE), India. One of the authors (R.K.M) thanks the University Grants Commission (UGC) for providing the Research Scholarship.

Figure Captions

Figure 1. The radial distribution function $`g(r)`$ of the system at $`T^{}=0.67`$, the lowest temperature investigated.

Figure 2. Temperature dependence of the shear viscosity ($`\eta `$). (a) ln($`\eta `$) as a function of inverse of temperature ($`1/T^{}`$). The simulating values given by the solid circles show super-Arrhenius behavior. The dashed line gives a guideline to the Arrhenius behavior. (b) ln($`\eta `$) against 1/($`TT_o^\eta `$). The solid circles are again represents the simulation results and the VFT fit is shown by the solid line. $`T_o^\eta `$ is found to be 0.57. The slope ($`E_\eta `$) and intercept (ln $`A_\eta `$) obtained from the fit are 0.81 and 0.69, respectively. For details, see the text.

Figure 3. Temperature dependence of the diffusion coefficients for the smallest ($`D_1`$) and largest ($`D_{10}`$) particles. (a) ln($`D_1`$) is plotted against 1/($`TT_o^{D_1}`$). (b) ln($`D_{10}`$) against 1/($`TT_o^{D_{10}}`$). Solid circles are the simulation results and the solid lines are the VFT fit. The critical temperatures $`T_o^{D_1}`$ and $`T_o^{D_{10}}`$ obtained from the VFT fit are 0.46 and 0.50, respectively.

Figure 4. The critical temperature ($`T_o^{D_i}`$) obtained from the VFT fit to the different subsets of particles as a function of the mean radius ($`\overline{R}_i`$) of the subsets (in units of $`\overline{\sigma }`$). Note that $`T_o^{D_i}`$ increases with the size of the particles.

Figure 5. The diffusion coefficients ($`D_i`$) as a function of 1/$`\overline{R}_i`$ at $`T^{}=0.67`$. The dashed line represents the Stokes-Einstein relation (Eq. 2) with the stick boundary condition $`C=6`$. The viscosity ($`\eta `$) value is taken from the present simulations. Solid circles are the simulated values and the solid line is the cubic polynomial fit in 1/$`\overline{R}_i`$. The fit parameters are as follows: $`D_i=0.00110.00132(1/\overline{R}_i)+0.000442(1/\overline{R}_i)^2`$. It clearly shows a highly nonlinear size dependence and a marked deviation from the Stokesian behavior of the diffusion.

Figure 6. The van Hove self-correlation function $`G_s(r,t)`$ as a function of the particle displacements $`r`$ (in units of $`\overline{\sigma }`$) at $`T^{}=0.67`$ for different values of time t (in units of $`\tau =\sqrt{\overline{m}\overline{\sigma }^2/ϵ}`$ = 2.2 ps for argon units). (a) For the smallest size particles (subset 1). The occurrence of the second peak at $`r1.0\overline{\sigma }`$ indicates the single particle hopping. (b) For the largest size particles (subset 10). Here also a second peak corresponds to single particle hopping develops but at relatively longer times.

Figure 7. The self-intermediate scattering function $`F_s(k,t)`$ for $`T^{}=0.67`$ is shown with a shift in the time origin to $`t_o=1.0`$, and normalized to the value at $`t_o`$, for a fixed value of $`k\overline{\sigma }=2\pi `$. This transformation is a convenient way to eliminate the Gaussian dependence at short time. (a) For the smallest size particles (subset 1). (b) For the largest size particles (subset 10). Open circles represents the simulation results and the solid lines are the stretched exponential fit (Eq. 6). The time constants ($`\tau _1`$ and $`\tau _{10}`$) and the exponents ($`\beta _1`$ and $`\beta _{10}`$) obtained from the fits are $`\tau _1242`$, $`\beta _10.49`$, $`\tau _{10}717`$, and $`\beta _{10}0.64`$.

Figure 8. (a) Projections into x-y plane of the trajectory of a typical smallest size particle over a time interval $`t=500\tau `$. (b) Projections into x-y plane of the trajectory of a typical largest size particle over a time interval $`t=2000\tau `$. Note that the time (t) is scaled by $`\tau =\sqrt{\overline{m}\overline{\sigma }^2/ϵ}`$ it is 2.2 ps if argon units are assumed. For detailed discussion, see the text.