text:
# Storing and releasing light in a gas of moving atoms

## Abstract

We propose a scheme of storing and releasing pulses or cw beams of light in a moving atomic medium illuminated by two stationary and spatially separated control lasers. The method is based on electromagnetically induced transparency (EIT) but in contrast to previous schemes, storage and retrieval of the probe pulse can be achieved at different locations and without switching off the control laser.

Over the past couple of years there has been a great deal of interest in storing and releasing light pulses in atomic gases exhibiting electromagnetically induced transparency (EIT) Fleisch:00 ; Hau:01 ; Lukin:01 ; Scully:01 ; Juz:02 ; Fleisch:02 ; Scully:02 . EIT is a phenomenon in which a weak probe pulse travels slowly and almost without dissipation in a resonant medium controlled by another laser Harris:90 ; Harris:91 ; Arimondo:96 ; Harris:97 ; Scully:book ; Bergman:98 ; Hau:99 . Following the proposal of Ref.Fleisch:00 , storage and release of a probe pulse has been demonstrated Hau:01 ; Lukin:01 ; Hemmer:02  by dynamically changing the intensity of the control laser. In the read-out process it is possible to alter certain properties of the stored probe pulse, such as its carrier frequency, propagation direction or pulse shape Juz:02 ; Scully:02 .

Recently an application of the light-storage scheme to generate continuous beams of atoms in nonclassical quantum states was proposed Gong . In that scheme the probe field propagates co-parallel to a beam of atoms in a spatially varying control field. This corresponds in the rest frame of the atoms to an explicitly time dependent control laser and thus allows for a complete and loss-free adiabatic transfer of probe-field excitations to the matter wave Fleisch:02 . We consider here another situation in which both the control and the incoming probe beams are perpendicular to the atomic motion. Such a scheme is better suited to eliminate effects of Doppler broadening and thus to alleviate limitations to the allowed velocity spread of the atomic beam. Furthermore in contrast to previous setups Fleisch:00 ; Hau:01 ; Lukin:01 ; Scully:01 ; Juz:02 ; Fleisch:02 ; Scully:02 , storing and releasing of a probe beam can now be achieved with a pair of stationary control lasers, i.e. there is no need to switch ”off” and ”on” a control laser at precise times. This is advantageous if one does not know the exact time of arrival of the probe pulse. Thus the present set-up could be used for example to store entangled photons generated by a continuous optical parametric oscillator below threshold and frequency selected by a cavity Schori02 ; Lvovsky01 . Finally storage and release of the probe field are spatially separated in the present scheme.

Consider a stream of atoms moving along the $`x`$ axes as shown in fig.1 (bottom part). The atoms are characterized by two hyper-fine ground levels $`g`$ and $`q`$, as well as an electronic excited level $`e`$ (top part of fig.1). The atoms in different internal states are described by the field-operators $`\mathrm{\Psi }_g\mathrm{\Psi }_g(𝐫,t)`$, $`\mathrm{\Psi }_q\mathrm{\Psi }_q(𝐫,t)`$ and $`\mathrm{\Psi }_e\mathrm{\Psi }_e(𝐫,t)`$ obeying Bose-Einstein or Fermi-Dirac commutation relations (depending on the type of atoms). The matter fields interact with several light fields propagating along the $`z`$-axis. A strong classical control laser centered at $`x=x_1`$ drives the transition $`|e|q`$, whereas a weaker quantum probe field (also entering the medium at $`x=x_1`$) is coupled with the transition $`|g|e`$. Finally there is a second spatially separated control laser (centered at $`x=x_2`$) that is used to release the stored probe pulse.

The first control laser has a frequency $`\omega _c`$, a wave vector $`𝐤_c=\widehat{𝐳}k_c`$, and a Rabi frequency $`\mathrm{\Omega }(𝐫,t)=\mathrm{\Omega }_1\left(x\right)e^{i\left(\omega _ctk_cz\right)},`$ where the $`y`$ dependence of the amplitude $`\mathrm{\Omega }_1\left(x\right)\mathrm{\Omega }_1`$ will be kept implicit. The same applies to other parameters, such as the mean atomic number density $`nn\left(z\right)`$.

The probe beam is described by the electric field operator: $`𝐄(𝐫,t)=\widehat{𝐞}\left(\mathrm{}ck/2\epsilon _0\right)^{1/2}(𝐫,t)e^{i\left(\omega tkz\right)}+H.c.`$ where $`(𝐫,t)`$ is the slowly varying amplitude, $`\omega =ck`$ is the central frequency of probe photons, $`𝐤=\widehat{𝐳}k`$ is the wave vector and $`\widehat{𝐞}\widehat{𝐳}`$ is the unit polarization vector. The dimension of $``$ is such that the operator $`^{}`$ represents the number density of probe photons.

The atoms (initially in the ground level $`g`$) are moving along the $`x`$ axes with an average speed $`𝐯_0=\mathrm{}𝐤_0/m`$ and the kinetic energy $`\omega _{at}=\mathrm{}^2k_0^2/2m`$. The atomic velocities are assumed to be spread over a narrow range around $`v_0`$, so that we can introduce slowly varying atomic amplitudes as: $`\mathrm{\Phi }_g=\mathrm{\Psi }_ge^{i\left(\omega _{at}t𝐤_0𝐫\right)}`$, $`\mathrm{\Phi }_e=\mathrm{\Psi }_ee^{i\left(\omega _{at}+\omega \right)ti\left(𝐤_0+𝐤\right)𝐫}`$ and $`\mathrm{\Phi }_q=\mathrm{\Psi }_qe^{i\left(\omega _{at}+\omega \omega _c\right)ti\left(𝐤_0+𝐤𝐤_c\right)𝐫}`$.

Consider first storing of the probe by the first control laser. The following equations hold for the slowly varying electromagnetic and matter field operators:

$`\left({\displaystyle \frac{}{t}}+c{\displaystyle \frac{}{z}}\right)`$ $`=`$ $`ig\mathrm{\Phi }_g^{}\mathrm{\Phi }_e,`$ (1)
$`\left({\displaystyle \frac{}{t}}+v_0{\displaystyle \frac{}{x}}\right)\mathrm{\Phi }_g`$ $`=`$ $`ig^{}\mathrm{\Phi }_e,`$ (2)
$`\left({\displaystyle \frac{}{t}}+i\mathrm{\Delta }+v_0{\displaystyle \frac{}{x}}+v_r{\displaystyle \frac{}{z}}\right)\mathrm{\Phi }_e`$ $`=`$ $`i\mathrm{\Omega }_1\mathrm{\Phi }_q+ig\mathrm{\Phi }_g,`$ (3)
$`\left({\displaystyle \frac{}{t}}+i\delta +v_0{\displaystyle \frac{}{x}}+\mathrm{\Delta }v_r{\displaystyle \frac{}{z}}\right)\mathrm{\Phi }_q`$ $`=`$ $`i\mathrm{\Omega }_1\mathrm{\Phi }_e,`$ (4)
where $`g=\mu \sqrt{ck/2\epsilon _0\mathrm{}}`$ characterizes the strength of the radiation-matter coupling, and $`\delta =\omega _q+\omega _c\omega \omega _g+\mathrm{\Delta }\omega _r`$ and $`\mathrm{\Delta }=\omega _e\omega \omega _g+\omega _r,`$ are the two- and single-photon detunings, $`v_r=\mathrm{}k/m`$ is the atomic recoil velocity due to absorption of a probe photon, $`\mathrm{}\omega _r=\mathrm{}k^2/2m`$ is the recoil frequency and $`\mathrm{\Delta }\omega _r=\mathrm{}\left(kk_c\right)^2/2m`$. The probe and control beams are assumed to be copropagating. In this case the overall recoil velocity $`\mathrm{\Delta }v_r=\mathrm{}\left(kk_c\right)/m`$ is small and can be neglected in eq.(4). Dissipation of the excited state $`|e`$ can be included in (3) replacing $`\mathrm{\Delta }`$ by $`\mathrm{\Delta }i\gamma `$ and adding the appropriate noise operator.

Assume that the probe field $``$ is sufficiently weak so that one can disregard the depletion of the gound level $`|g`$. Neglecting the last term in eq.(2), one has

$$\frac{}{t}\mathrm{\Phi }_g(𝐫,t)=v_0\frac{}{x}\mathrm{\Phi }_g(𝐫,t).$$
(5)
It is convenient to introduce the field operators describing annihilation of atomic and spin excitations (excitons):

$$\psi _u\psi _u(𝐫,t)=n^{1/2}\mathrm{\Phi }_g^{}\mathrm{\Phi }_u$$
(6)
with $`u=e,q`$, where $`nn(z)=\mathrm{\Phi }_g^{}\mathrm{\Phi }_g`$ describes the spatial profile of the number-density of ground-state atoms. Exploiting relation (5), the equations of motion for the new operators have the form of eqs.(1), (3) and (4) subject to the replacement $`\mathrm{\Phi }_e\psi _e`$, $`\mathrm{\Phi }_q\psi _e`$, and $`\mathrm{\Phi }_gn^{1/2}`$. Note that the exciton operators $`\psi _e`$ and $`\psi _q`$ obey approximately Bose-commutation relations even through the constituent atoms may be fermions.

Let us consider the case of exact two-photon resonance ($`\delta =0`$). This is well justified as the two-photon Doppler shift due to longitudinal and transversal velocity spreads of the atoms is strongly diminished due to the chosen geometry. Neglecting terms containing $`\psi _e`$ and $`\dot{\psi }_e`$ in eq. (3), one arrives at the adiabatic approximation relating $`\psi _q`$ to the electric field amplitude $`(𝐫,t)`$ as:

$$\psi _q(𝐫,t)=gn^{1/2}\stackrel{~}{}(𝐫,t)$$
(7)
with $`\stackrel{~}{}(𝐫,t)=(𝐫,t)/\mathrm{\Omega }_1(x)=\psi _q(𝐫,t)/gn^{1/2}`$ being an auxiliary field. Assuming a stationary flow of atoms in the $`x`$-direction, the atomic density $`nn\left(z\right)`$ does not depend on $`x`$, and thus eqs. (4) and (7) yield

$$\psi _e(𝐫,t)=i\frac{gn^{1/2}}{\mathrm{\Omega }_1}\left(v_0\frac{}{x}+\frac{}{t}\right)\stackrel{~}{}(𝐫,t).$$
(8)
On the other hand, the Rabi frequency $`\mathrm{\Omega }_1\mathrm{\Omega }_1\left(x\right)`$ is $`z`$ independent, so eqs. (1) and (8) lead to the following equation for $`\stackrel{~}{}(𝐫,t)`$:

$$\left(\frac{}{t}+\frac{c}{1+n_g}\frac{}{z}+\frac{v_0n_g}{1+n_g}\frac{}{x}\right)\stackrel{~}{}(𝐫,t)=0,$$
(9)
with $`n_g(x,z)=g^2n(z)/\mathrm{\Omega }_1^2(x)`$ being the group index. Assuming slow propagation, the group index is much larger than unity, so that eq. (9) simplifies to

$$\left(\frac{}{t}+\stackrel{~}{v}_g\left(z\right)a\left(x\right)\frac{}{z}+v_0\frac{}{x}\right)\stackrel{~}{}(𝐫,t)=0,$$
(10)
where $`v_g\stackrel{~}{v}_g\left(z\right)a\left(x\right)=c/n_g(x,z)`$ is the group velocity ($`v_gc`$). The dimensionless quantity $`a\left(x\right)a_1\left(x\right)=\left[\mathrm{\Omega }_1\left(x\right)/\mathrm{\Omega }_1\left(x_1\right)\right]^2`$ characterizes the spatial shape of the first control laser centered at $`x=x_1=0`$. The $`z`$ dependence of $`\stackrel{~}{v}_g\left(z\right)v_g(z,x_1)`$ emerges through the atomic density $`nn\left(z\right)`$. It is noteworthy that $`(𝐫,t)=\psi _q(𝐫,t)\sqrt{v_g/c}`$, so the ratio between the number density of photons and the spin excitations is given by the relative group velocity $`v_g/c`$. In other words, the slowly propagating probe beam is made of the EIT (dark state) polaritons Fleisch:00 ; Juz:02 ; Fleisch:02  comprising predominantly the spin excitations.

Assume that the spatial width of the first control beam $`\mathrm{\Delta }x_{c1}`$ is much larger than that of the incomming probe beam $`\mathrm{\Delta }x_p`$ centered at $`x=x_1=0`$. Under this condition, the solution of eq.(10) can be expressed in terms of the incoming electric field at an entry point $`z=z_1`$ as

$$\stackrel{~}{}(𝐫,t)=\frac{[\xi (x,z),z_1,\tau (t,x,z)]}{\mathrm{\Omega }_1\left(x_1\right)},$$
(11)
where $`\xi (x,z)=_0^xa\left(x^{}\right)𝑑x^{}v_0_{z_1}^z\left[\stackrel{~}{v}_g\left(z^{}\right)\right]^1𝑑z^{}`$, and $`\tau (t,x,z)=t+\xi (x,z)/v_0x/v_0`$ obey the proper boundary conditions: $`\xi \left(x,z=z_1\right)x`$ and $`\tau \left(t,x,z=z_1\right)t`$ for the incoming probe field.

Using the definition of $`\stackrel{~}{}`$ and eq.(11), one finds the temporal and spatial behavior of the operators for electric field and spin excitations:

$`(𝐫,𝐭)`$ $`=`$ $`\sqrt{|a\left(x\right)|}[\xi (x,z),z_1,\tau (t,x,z)],`$ (12)
$`\psi _q(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{gn^{1/2}}{\mathrm{\Omega }_1\left(x_1\right)}}[\xi (x,z),z_1,\tau (t,x,z)].`$ (13)
Equations (12) and (13) define the electric and spin components of the EIT polariton. A set of trajectories for such a polariton in the $`xz`$ plane is given by $`\xi (x,z)=\xi _0`$, with $`\xi _0=0`$ corresponding to the central trajectory. Due to the motion of the atoms parallel to the $`x`$ axis, the polariton is dragged into that direction. Following the spatial profile $`a(x)`$ of the coupling beam, its velocity component in the $`z`$ direction, $`\stackrel{~}{v}_g(z)a(x)`$, is further reduced. As soon as $`\stackrel{~}{v}_g(z)a(x)`$ becomes less than $`v_0`$ the flow of excitation is predominatly determined by the velocity of the atoms $`v_0`$ in the $`x`$ direction. If the atomic beam is optically thick in the $`z`$direction ($`\xi (\mathrm{},z_{max})<0`$), the propagation direction of the polariton may be completely converted from $`\widehat{𝐳}`$ to $`\widehat{𝐱}`$. The probe field is then stored in the form of a spin excitation that moves with the atoms parallel to the $`x`$-axis and is centered at $`z=z_{\mathrm{}}`$, which is a solution of $`\xi \left(x\mathrm{},z_{\mathrm{}}\right)=0`$. This is illustrated in figs.2 and 3 where we have shown the amplitude of the electric field at different instances of time. Also shown is the effect of the second regeneration laser, which will be discussed later on.

As one can see from the figures, the dragging of the polariton along with the moving atoms leads to a deformation of the excitation. For large values of $`x`$, such that $`a(x)0`$, one has in the vicinity of $`z=z_{\mathrm{}}`$

$`\xi \left(x\mathrm{},z\right)\left(zz_{\mathrm{}}\right){\displaystyle \frac{v_0}{\stackrel{~}{v}_g(z_{\mathrm{}})}},`$ (14)
$`\tau (t,x\mathrm{},z){\displaystyle \frac{zz_{\mathrm{}}}{\stackrel{~}{v}_g(z_{\mathrm{}})}}{\displaystyle \frac{xx_0}{v_0}},`$ (15)
where $`x_0x_0\left(t\right)=v_0t`$ is the ”center of mass” of a stored polariton moving in the $`x`$ direction. Thus, if $`\mathrm{\Delta }x_p`$ and $`2\mathrm{\Delta }\tau _p`$ denote the half-width (in $`x`$ direction) and the duration of the input probe pulse, the dimensions of the spin wave are in the limit $`x\mathrm{}:`$

$`\mathrm{\Delta }x_sv_0\mathrm{\Delta }\tau _p,\mathrm{\Delta }z_s\mathrm{\Delta }x_p{\displaystyle \frac{\stackrel{~}{v}_g(z_{\mathrm{}})}{v_0}},`$ (16)
where we have assumed that $`v_0\mathrm{\Delta }x_p/\mathrm{\Delta }\tau _p`$. Since the spatial profile of the polariton in the $`y`$ direction is not changed, the transfer from a electromagnetic to a matter-wave pulse is associated with a spatial compression of the excitation volume by a factor $`\stackrel{~}{v}_g(z_{\mathrm{}})/c`$. From eq.(16) one can also easily obtain a necessary condition for a complete transfer of the electromagnetic excitation to the atomic beam: Since $`\mathrm{\Delta }z_s`$ should be less than the half width of atomic beam $`\mathrm{\Delta }z_{\mathrm{atom}}`$, one finds a minimum ratio of $`v_0`$ to the group velocity at the center

$$\frac{v_0}{\stackrel{~}{v}_g(z_{\mathrm{}})}\frac{\mathrm{\Delta }x_p}{\mathrm{\Delta }z_{\mathrm{atom}}}.$$
(17)
If this condition is fulfilled a pulse is completely stored within the atom beam and does not exit on the back side.

Consider now the regeneration of the probe beam by a spatially separated, i.e. nonoverlapping, second control laser centered at $`x=x_2`$ and characterized by a Rabi frequency $`\mathrm{\Omega }_2\left(x\right)`$. In this case, the previous equations (10) - (14) describe the propagation of a probe beam within the entire system, subject to the following replacement: $`a\left(x\right)a\left(x\right)=a_1\left(x\right)\pm a_2\left(x\right)`$ where $`a_2\left(x\right)=\left[\mathrm{\Omega }_2\left(x\right)/\mathrm{\Omega }_2\left(x_1\right)\right]^2`$ characterizes the shape of the second control laser. The upper (lower) sign corresponds to the case where the second control laser propagates in the same (opposite) direction, as compared to the first one footn:2 . This is because the radiative group velocity of the regenerated probe beam changes sign in the latter case. It is noteworthy that the reversed probe beam experiences a slight shift out of the EIT resonance Juz:02 . Yet such a shift can be neglected in the case where the initial control and probe beams are copropagating. Note also that the previously considered regeneration of a probe beam is due to the temporal switching of the control laser Fleisch:00 ; Hau:01 ; Lukin:01 ; Scully:01 ; Juz:02 ; Fleisch:02 ; Scully:02 , whereas now the regeneration is induced by a second spatially separated control laser that can be stationary.

Finally, let us analyze the output field at a spatial point where $`z=z_2`$. The maximum electric and spin fields are then concentrated at $`x=x_2`$, for which $`\xi (x_2,z_2)=0`$. In the vicinity of this point, one has $`\xi (x,z_2)\pm a_2\left(x_2\right)\left(xx_2\right)`$ i.e. the output field has a width $`\mathrm{\Delta }x_{p2}=\mathrm{\Delta }x_{p1}/a_2\left(x_2\right)`$.

In summary, we have investigated a different scheme of storing and releasing a beam of probe light in a moving atomic medium illuminated by two spatially separated control lasers depicted in fig.1. Beyond the area illuminated by the first control laser, the probe beam transforms into a beam of pure spin excitations moving along the $`x`$ axis, as one can see from figs.2,3. The regeneration of the probe beam is accomplished by applying the second continuous control laser. Depending on the direction of the latter, the restored probe beam moves either parallel (fig.2), or antiparallel (fig.3) to the initial probe beam.

In contrast to the previous schemes Fleisch:00 ; Hau:01 ; Lukin:01 ; Scully:01 ; Juz:02 ; Fleisch:02 ; Scully:02 , storing and releasing of a probe beam can be now accoplished without switching ”off” and ”on” of a control laser. This is advantageous if one does not know the exact time of arrival of the probe photons, e.g. if the latter are created via spontaneous processes.


format: text

text:
# Untitled Document

The alternative to the incompressible fractional charge in quantum Hall effect: Comments on Laughlin and Schrieffer’s papers. Keshav N. Shrivastava School of Physics, University of Hyderabad,
Hyderabad 500046, India

Laughlin has found “exactly” the wave function which is ascribed to an excitation of fractional charge, such as e/3. We find that the exactness of the wave function is not destroyed by changing the charge to some other quantity, such as the magnetic field. Thus e/3 and H can be replaced by e and H/3. Therefore, the wave function need not belong to a quasiparticle of charge e/3. Corresponding author: keshav@mailaps.org
Fax: +91-40-3010145.Phone: 3010811.

1. Introduction

An effort is made to see if there is an alternative to the interpretation of the fractional charge in the quantum Hall effect. Laughlin has written the wave function for the quasiparticle of fractional charge. We examine to see if the wave function can be interpreted to describe integer charge and the blame of the fraction can be thrown on some other quantity. Laughlin1 has proposed a variational ground state which is a new state of matter, a quantum fluid, the elementary excitations of which are fractionally charged. The correctness of the wave functions is verified by direct numerical diagonalization of the many-body Hamiltonian. We wish to examine the alternative interpretation of the wave function without destroying the exactness of the calculation. In particular, we examine whether some other variable can take the blame of the fraction instead of the charge. If we bring a flux quantum near a solenoid, according to Laughlin, charge is affected whereas according to the correct answer, the field is affected. The charge density should automatically give the correct electric and magnetic vectors so that the magnetic field should be automatically correct. There are two variables, the charge and the field in the flux quantization. Therefore, Laughlin’s throwing the blame only on the charge is not correct and the automatic correction to the electric and magnetic vectors does not occur.

We find that several quantities arise as a factor of the charge so that the effect of the fraction need not be thrown on the charge. Since the charge multiplies the magnetic field, the factor arising in the field can also be read as a factor of the charge, then there is no trouble in defining a fractional charge. However, if charge is corrected without correcting the magnetic field, apparently there is no trouble, except that the choice of the variale is not unique.

In this comment, we explain the various choices available without disturbing the exactness and point out that field is better than the charge.

2. Theory

Let us consider a two-dimensional electron gas in the $`xy`$ plane subjected to a magnetic field along z direction. The eigen states of the single-body Hamiltonian can be written as, $`|m,n>`$. The cyclotron frequency is, $`\mathrm{}\omega _c`$=$`\mathrm{}(eH_o/mc)`$. The magnetic length is $`a_o`$=($`\mathrm{}/m\omega _c`$)1/2 which upon substituting the cyclotron frequency becomes,

$$a_o=(\frac{\mathrm{}c}{eH_o})^{1/2}.$$
(1)
The energy levels of the type of a harmonic oscillator are produced when the applied field is along the z direction,

$$H|m,n>=(n+\frac{1}{2})|m,n>.$$
(2)
The wave function of the lowest Landau level is written as a function of $`z=x+iy`$ with $`|m>`$ as an eigen state of angular momentum with eigen value $`m`$. The many-body Hamiltonian consists of, the kinetic energy with the vector potential included in the momentum, a potential due to the positively charged nuclei, V(zj), and the Coulomb repulsion between electrons. The wave functions composed of states in the lowest Landau level which describe the angular momentum $`m`$ about the center of mass are of the form,

$$\psi =(z_1z_2)^m(z_1+z_2)^nexp[(\frac{1}{4})(|z_1|^2+|z_2|^2)].$$
(3)
Laughlin generalized this observation to $`N`$ particles by writing product of Jastrow functions,

$$\psi =\{\mathrm{\Pi }_{j<k}f(z_jz_k)\}exp(\frac{1}{4}\mathrm{\Sigma }_l|z_l|^2)$$
(4)
which minimizes the energy with respect to $`f`$. If $`\psi `$ is antisymmetric $`f(z)`$ must be an odd function, $`f(z)`$ = zm with m = odd. To determine which $`m`$ minimizes the energy, we write,

$$|\psi _m|^2=exp(\beta \varphi )$$
(5)
where $`\beta `$=1/m and $`\varphi `$ is the classical potential energy which describes a system of $`N`$ identical particles of charge $`Q`$=m with the neutralizing back ground charge density, $`\sigma `$=(2$`\pi `$ a$`{}_{}{}^{2}{}_{o}{}^{}`$)-1 per unit area. This is the classical one-component plasma (OCP). The solution of which is well known. For $`\mathrm{\Gamma }`$=2$`\beta `$Q2=2m $`>`$ 140, the OCP is a hexagonal crystal and fluid otherwise. Laughlin’s wave function describes a fluid of density,

$$\sigma _m=\frac{1}{m(2\pi a_o^2)}$$
(6)
which minimizes the energy. The charge density generated by (6) should be equal to that of the background charge so that there is overall charge neutrality.

If we define a new value of $`a_o`$, then $`ma_o^2`$ is replaced by $`a_{new}^2`$. Then the effect of $`m`$ is not on the charge but it is in the distance. That means that the effect of $`m`$ can occur on, (i)h, (ii)c,(iii) e or (iv) Ho. Thus there are at least four candidates to absorb the effect of $`m`$. Let us eliminate the Planck’s constant and the velocity of light. Then it is possible to affect either $`e`$ or $`H_o`$. It will surely be interesting if the effect of $`m`$ is thrown on the velocity of light. Then,

$$a_{new}^2=\frac{\mathrm{}(mc)}{eH_o}$$
(7)
where for $`m`$=3, the quasiparticles will travel with the velocity of 3$`c`$ which is faster than light. In the expression (1), the velocity of light $`c`$ is the value in vacuum so that the particles with velocity 3$`c`$ will be travelling faster than light which will be noncausal or causality violating. Another interpretation is obtained by throwing the value of $`m`$ on the magnetic field,

$$a_{new}^2=\frac{\mathrm{}c}{e(H_o/m)}.$$
(8)
Therefore, the number $`m`$ of Laughlin can be absorbed as a factor of magnetic field such as,

$$gH_o=H_o/m.$$
(9)
It is perfectly allowed to devide the magnetic field by $`m`$ to define a new field. Out of the four options available, i.e., h, c, e and Ho, Laughlin selected $`e`$ so that, the effective charge becomes,

$$e_{eff}=e/m$$
(10)
where h, c and Ho are kept constant and $`m`$ is an odd integer. Therefore, Laughlin’s choice of variables is not unique. Let us ask, if some one had selected the velocity of light, out of the four variables, then what would have happened? The answer to this question is that we would have obtained the particles faster than light.

3. Exactness

The projection of $`\psi _m`$ for three and four particles onto the lowest energy eigen state of angular momentum 3$`m`$ for m=3 and m=5 states is about 0.98 which shows the exactness of the wave functions. Therefore, the exactness of the calculation is not affected by changing the variable from $`e`$ to $`H`$, h or c. The calculation is exact for any of the four variables. The total energy per particle can be written in terms of the distribution function of the one-component plasma (OCP),

$$U_{tot}=\pi _o^{\mathrm{}}\frac{e^2}{r}[g(r)1]r𝑑r(4/3\pi 1)2e^2/R$$
(11)
where integration domain is a disk of radius R=($`\pi \sigma _m`$)-1/2. At $`\mathrm{\Gamma }`$=2, g(r) = 1- exp\[-(r/R)2\] which shows that for $`m<`$9, the Utot is deeper than for charge density waves. The excitations of $`\psi _m`$ are created by piercing the fluid at $`z_o`$ with an infinitely thin solenoid and passing through it a flux quantum $`\mathrm{\Delta }\phi `$=hc/e, adiabatically. This treatment does not treat the magnetic field correctly. Here again,

$$field\times area=hc/e$$
(12)
so that the variable is e$`\times `$field$`\times `$area and hence the entire blame can not be thrown on $`e`$. It may appear that changing the $`e`$ will automatically take care of the electric vector E and the magnetic field H of the electromagnetic field but there is H in the flux quantization which is unchanged. A change in the charge density should change both E and H. According to the adiabatic approximation, the flux quantum produces the quasiparticles. If this approximation is not valid, field will be produced not the quasiparticles. In fact the variables, $`h`$ and $`c`$ are also very good candidates. Increased value of the Planck’s constant will increase the cyclotron absorption energy and increased $`c`$ will give tycheons: noncausal faster than light particles. The effect of passing a flux quantum is to change the single-body wave function from,

$$(zz_o)^mexp(\frac{1}{4}|z|^2)to(zz_o)^{m+1}exp(\frac{1}{4}|z|^2).$$
(13)
An approximate representation of these states is chosen by Laughlin for the quasielectron and the quasihole. Laughlin writes $`|\psi ^{+z_o}|^2`$ as exp(-$`\beta \varphi `$) with $`\beta `$=1/$`m`$ and $`\varphi `$=$`\varphi `$-2$`\mathrm{\Sigma }_lln|z_lz_o|`$ where $`\varphi `$ describes an OCP interacting with a phantom point charge at $`z_o`$. The plasma screens this phantom by accumulating an equal and opposite charge near $`z_o`$. The particle charge in the plasma is 1, rather than $`m`$ so the accumulated charge is 1/m, i.e., meeff = 1. The energy required to create a particle of Debye length $`a_o/\sqrt{2}`$ is given by,

$$\mathrm{\Delta }_D=\frac{\pi }{4\sqrt{2}}\frac{e^2}{m^2a_o}.$$
(14)
Note that the quantity which enters is again $`m^2a_o`$ and not just $`e/m`$ so that the charge can be changed from $`e`$ to $`e/m`$ or the charge may be kept constant at $`e`$ and $`a_o`$ changed to $`m^2a_o`$. Therefore, two equivalent possibilities exist. One is to define the effective charge, eeff=$`e/m`$ and the other is to keep $`e`$ unchanged and change $`a_o`$ to $`a_{eff}`$= $`m^2a_o`$. The state described by $`\psi _m`$ is incompressible because compressing it is equal to injecting particles and particles carry charge so the fractional charge will be destroyed. The incompressibility is not built in the hamiltonian and is extranuously imposed. It is easy to impose an external boundary condition that there is incompressibility but in the expression (14) such an incompressibility can not be enforced. Therefore, the incompressibility enforcement is not contained in the theory. The compressibility creates sound waves so the incompressibility causes the sound to be absent. This is equivalent to eliminating phonons, which has not been done, so the resistivity can touch zero value. In the BCS theory of superconductivity, the zero resistance is obtained not by introducing incompressibility but by eliminating phonons which makes the electrons attactive. We do not discuss the impurities or rotations because these are not contained in our hamiltonian. The Hall conductance is (1/m)e2/h so that the effective charge is eeff=(1/m)e and the factor e/h is caused by the units. The origin of this effective charge is not clear because of the other candidates and for m=3, the charge e/3 is not generated by the hamiltonian discussed by Laughlin but the antisymmetry surely requires that m is odd, such as 1, 3 , 5, 7, … etc. Therefore, odd aspect is a result of antisymmetry. Laughlin does not use spin so that there is spin-charge decoupling automatically built in the calculation. The $`\mathrm{𝐬𝐩𝐢𝐧𝐥𝐞𝐬𝐬}`$ electrons may give unphysical results. Below eq.(13) it is argued that introducing a flux quantum $`\mathrm{\Delta }\varphi `$=hc/e results into accumulation of charge 1/$`m`$. In fact, the charge need not accumulate and only the field is modified as in (8). Therefore, the factor 1/$`m`$ does not determine the charge of the excitations. In a different context Laughlin has pointed out that the bulk modulus of the Hartree-Fock ground state was $`\mathrm{𝑧𝑒𝑟𝑜}`$. The bulk modulus is actually finite, as is usually the case for Jastrow-type trial wave functions for helium as admitted in an errata2.

4. Missing charge.

If we agree to the choice of the charge and change it from $`e`$ to $`e`$/3, then what happened to the remaining charge of 2/3? If we say that $`e`$ remains $`e`$ and $`H`$ changes to $`H`$/3, then we need not look for the missing charge. The wave function (4) does not conserve charge. Laughlin does not provide any prescription to find the missing charge.

5. Conclusions

The exact wave function, the excitations of which are fractionally charged does not determine the quasiparticle charge uniquely. The blame of the fraction need not be thrown on the quasiparticle charge. The introduction of the alternative choice of variable, instead of charge, does not destroy the exactness of the calculation. The magnetic field of the flux quantization condition has been left out.

It is clear that Laughlin’s wave function can not explain the experimental observations of the quantum Hall effect. It may be pointed out that the composite fermion model (CF) which requires that the flux quanta be attached to the electron is also not correct.3-5

In the expression (13) of Arovas et al6 the quantity $`e^{}A_\varphi `$ can be replaced by (e/m)Aϕ. Now, multiply $`A_\varphi `$ by $`g`$ so that the quantity of interest becomes (e/m)gAϕ with g=1. It is clear that (e/m) and (g/m) are not resolvable because (e/m)gAϕ is exactly equal to e(g/m)Aϕ. Therefore, whether the charge should be fractionalized or the unit flux should be corrected are not resolved. If charge is corrected, then fractional statistics will come otherwise the unit flux will be changed. This means that the fermion need not obey the “fractional statistics” and $`\varphi _o`$ will be modified.

The correct theory of the quantum Hall effect is given in ref.7.

6. References

1. R. B. Laughlin, Phys. Rev. Lett. 50, 1395(1983).
2. R. B. Laughlin, Phys. Rev. Lett. 60, 2677 (1988);61, 379(E)(1988).
3. M. I. Dyakonov, cond-mat/0209206.
4. B. Farid, cond-mat/0003064.
5. K. N. Shrivastava, cond-mat/0209666.
6. D. Arovas, J. R. Schrieffer and F. Wilczek, Phys. Rev. Lett. 53, 722 (1984).
7. K. N. Shrivastava, Introduction to quantum Hall effect,
Nova Science Pub. Inc., N. Y. (2002).

Note: Ref.7 is available from:
Nova Science Publishers, Inc.,
400 Oser Avenue, Suite 1600,
Hauppauge, N. Y.. 11788-3619,
Tel.(631)-231-7269, Fax: (631)-231-8175,
ISBN 1-59033-419-1 US$`\$69`$.
E-mail: novascience@Earthlink.net


text:
# Spin-polarized stable phases of the 2-D electron fluid at finite temperatures.

## Abstract

The Helmholtz free energy $`F`$ of the interacting 2-D electron fluid is calculated nonperturbatively using a mapping of the quantum fluid to a classical Coulomb fluid $`[`$Phys. Rev. Letters, 87, 206 (2002)$`]`$. For density parameters $`r_s`$ such that $`25<r_s`$, the fluid is unpolarized at all temperatures $`t=T/E_F`$, where $`E_F`$ is the Fermi energy. For lower densities, the system becomes fully spin polarized for $`t`$ smaller than $`0.35`$, and partially polarized for $`0.35<t<2`$, depending on the density. At $`r_s2530`$, and $`t0.35`$, an ”ambispin” phase where $`F`$ is almost independent of the spin polarization is found. These results support recent claims, based on quantum Monte Carlo results, for a stable, fully spin-polarized phase at $`T`$ = 0 for $`r_s`$ larger than $`2526`$.

Introduction. The physics of the uniform two-dimensional electron gas (2DEG) depends crucially on the “coupling parameter” $`\mathrm{\Gamma }`$ = (potential energy)/(kinetic energy) arising from the Coulomb interactions. The $`\mathrm{\Gamma }`$ for the 2DEG at $`T=0`$ and mean density $`n`$ is equal to the mean-disk radius $`r_s=(\pi n)^{1/2}`$ per electron. In the absence of a magnetic field or disorder effects, the density parameter $`r_s`$, the spin polarization $`\zeta `$ and the temperature $`T`$ are the only relevant physical variables. As the coupling constant increases, the 2DEG becomes more like a “liquid”, but here we continue to loosely call it an electron “gas”. A knowledge of the phase diagram at arbitrary $`r_s`$, $`\zeta `$ and $`T`$ is of crucial importance to a clear understanding of many intriguing phenomena (metal-insulator transitionrevmod , 2-D electrons in CuO layers, spin-dependent processes etc.) associated with the 2DEG.

A number of recent studiesvarsen ; atta ; prl2  have suggested that the 2DEG at $`T`$ =0 becomes fully spin-polarized for $`r_s>r_c`$ with $`r_c2530`$. However the energy difference between the unpolarized and polarized phases is very close to the statistical uncertainty of the methods used. We have completed a study of this transition at finite temperatures where the energy differences for a larger data base become more significant and reveal a rich phase diagram which is fully consistent with the claimed $`T=0`$ polarized-fluid regime. It turns out that the partially degenerate strongly correlated 2DEG supports partially spin-polarized phases. This is particularly interesting since the 3DEG has also been recently shown to have partially spin-polarized stable phasesichi3dphase

The finite-$`T`$ phase diagram is examined using a computationally simple, conceptually novel method for calculating the quantum-PDFs and other properties (e.g, static response). The method has been applied to the 2DEGprl2 , the 3DEGprl1 ; prb00 , and hydrogen fluidshyd , at arbitrary coupling, spin polarization and temperature. It is based on identifying a “quantum temperature” $`T_q`$ such that the correlation energy of the corresponding classical Coulomb fluid at $`T_q`$ is equal to that of the quantum fluid at $`T=0`$. After that, hypernetted-chain (HNC) methods are used for the classical fluid. This classical mapping of quantum fluids within the HNC was named the CHNC. The application of the CHNC method to the 2DEG required the inclusion of short-range correlation effects (“bridge terms”) going beyond the usual HNC approximation. The resulting CHNC procedure for the 2DEG proved to be remarkably accurateprl2 . In fact even an application of CHNC at $`T`$ = 0, and without bridge corrections was examined by Bulutay and Tanatar and found to be very usefulbalutay .

We use Hartree atomic units, and consider a unit area $`A`$ with $`n`$ particles per unit area, with $`n`$ = $`1/(\pi r_s^2)`$. The number density of each spin species is $`n_i=n(1+\sigma _i\zeta )/2`$, where i=1,2 would indicate up or, down spins with $`\sigma _i=\pm 1`$. The spin polarization $`\zeta `$ = $`(n_1n_2)/n`$. The total Helmholtz free energy per unit area, $`F`$, has the form $`F_0+F_x+F_c`$, where $`F_0`$, $`F_x`$, $`F_c`$ are the non-interacting, first-order exchange, and correlation, contributions to $`F`$. Thus all effects beyond first-order exchange are lumped into $`F_c`$, a convention used in density functional theory (DFT). The DFT exchange-correlation free energy would hence be $`F_{xc}`$ = $`F_x+F_c`$. In indicating spin-components we use the notation $`F_i^0,F_i^x`$, and $`F_i^c`$ etc. The canonical ensemble is used so that the Helmholtz free energy is directly determined, unlike in the approach where the grand-potential $`\mathrm{\Omega }`$ is determined in terms of the chemical potential which has to be inverted to obtain the density. This means $`F_0`$ and $`F_x`$ (being the first-order exchange) are determined in terms of the non-interacting chemical potential $`\mu _0`$, while $`F_c`$ has to be determined via a coupling constant integration, as discussed by us in ref prb00 ; pdw84  and in Dundrea et aldundrea .

The non-interacting free energy. Isihara and others have considered the expansion of the free energy about $`T=0`$ and obtained the leading terms in $`r_s`$, and $`T`$ isihara . However, the thermodynamic functions near $`T=0`$ have a particularly unfavourable structure for constructing such expansions. Since our calculations are nonperturbative, we evaluate $`F_0`$ etc., numerically, without expansion.

Using the the relation, $`F_0=E_0TS_0`$, where $`E_0`$ and $`S_0`$ are the noninteracting internal energy and the entropy, $`F_0`$ can be calculated. An independent evaluation is obtained from the relation $`F_0=n\mu _0E_0`$. Also, at $`T=0`$, we have, in Hartrees

$$F_i^0/n=E_i^0/n=(1+\sigma _i\zeta )^2/4r_s^2$$
(1)
Defining the Fermi function $`f(e)=1/(1+exp(\beta e\eta ))`$ where $`\beta =1/T`$, $`\eta =\mu ^0/T`$,

$`n_i`$ $`=`$ $`(T/2\pi )I_0(\eta _i)`$ (2)
$`E_i`$ $`=`$ $`(T^2/pi)I_1(\eta _i)`$ (3)
$`F_i^0`$ $`=`$ $`n_i\mu _i^0E_i^0`$ (4)
where the Fermi integral $`I_\nu (z)`$ is defined as

$$I_\nu (z)=_0^{\mathrm{}}\frac{dxx^\nu }{1+e^{xz}}$$
(5)
The $`I_0`$ integral can be evaluated explicitly to give:

$`\eta _i`$ $`=`$ $`log[exp(1/t_i)1],t_i=E_i^F/T`$ (6)
$`E_i^F`$ $`=`$ $`E_F(1+\sigma _i\zeta ),E_F=\pi n`$ (7)
The $`I_1`$ integral can be evaluated directly or reexpressed in terms of the dilogarithm function. We have evaluated $`F_0`$ by several methods; the dilog function was evaluated using the CERN library routine.

The exchange contribution to the free energy. The first-order unscreened exchange free energy $`F_x`$ consists of $`F_x^i`$, where $`i`$ denotes the spin species. At $`T=0`$ these reduce to the exchange energies:

$$E_i^x/n=\frac{8}{3\sqrt{\pi }}n_i^{1/2}$$
(9)
Hence,

$$E_x/n=(E_1^x+E_2^x)/n=\frac{8}{3\pi r_s}[x_1^{3/2}+x_2^{3/2}]$$
(10)
where $`x_1`$ and $`x_2`$ are the fractional compositions $`(1\pm \zeta )/2`$ of the two spin species. At finite $`T`$, we have:

$$F_i^x/E_i^x=\frac{3}{16}t_i^{3/2}_{\mathrm{}}^{\eta _i}\frac{I_{1/2}^2(u)du}{(\eta _iu)^{1/2}}$$
(11)
Here $`I_{1/2}`$ is the Fermi integral. The total exchange free energy is $`F_x=\mathrm{\Sigma }F_i^x`$. The accurate numerical evaluation of Eq. 11 requires the removal of the square-root singularity by adding and subtracting, e.g, $`I^2(|\eta |)/(v|eta|)^{1/2}`$ for the case where $`\eta `$ is negative, and $`v=u`$, and so on.

A real-space formulation of $`F_x`$ = $`F_1^x+F_2^x`$ using the zeroth-order PDFs fits naturally with the CHNC approach. Thus

$$F_x/n=n\frac{2\pi rdr}{r}\underset{i<j}{}h_{ij}^0(r)$$
(12)
Here $`h_{ij}^0(r)=g_{ij}^0(r)1`$. In the non-interacting system at temperature $`T`$, the antiparallel $`h_{12}^0`$, viz., $`g_{12}^0(r,T)1`$, is zero while

$$h_{11}^0(𝐫)=\frac{1}{n_i^2}\mathrm{\Sigma }_{𝐤_1,𝐤_2}n(k_1)n(k_2)e^{i(𝐤_1𝐤_2)𝐫}=[f(r)]^2$$

Here k,r are 2-D vectors and $`n(k)`$ is the Fermi occupation number at the temperature $`T`$. At $`T=0`$ $`f(r)=2J_1(k_ir)/k_r`$ where $`J_1(x)`$ is a Bessel function. We have numerically evaluated the exchange free energy by both methods,i.e, Eq. 11 and Eq. 12, as a numerical check.

Correlation free energy. The correlation contribution, $`F_c`$ is evaluated together with the exchange contribution $`F_x`$ via a coupling constant integration over the distribution functions as follows:

$$F_{xc}/n=_0^1𝑑\lambda n\frac{2\pi rdr}{r}\underset{i<j}{}h_{ij}(r,\lambda )$$
(13)
Thus we need the interacting PDFs $`g_{ij}(r,\lambda )`$ for each $`r_s`$, $`\zeta `$ and $`T`$ for a number of values (we have used 7 values) of the coupling constant $`\lambda `$. The PDFs were evaluated using the CHNC method, as described in Ref.\[prl2, \]. We present a brief summary of the CHNC method in the context of the present study.

Brief summary of the CHNC method. The essence of the CHNC method is to start from the quantum mechanical $`g^0(r)`$ of the non-interacting problem and build up the interacting $`g(r)`$ by classical methods.

The fluid of mean density $`n`$ contains two spin species with concentrations $`x_i`$ = $`n_i/n`$. While $`T`$ is the physical temperature of the 2DEG, we consider a classical fluid of temperature $`T_{cf}`$ = $`1/\beta _{cf}`$. Since the leading dependence of the energy on temperature is quadratic, we assume that $`T_{cf}`$= $`\sqrt{(T^2+T_q^2)}`$ where $`T_q`$ is the quantum correction. This is clearly valid for $`T=0`$ and for high $`T`$, and was justified in more detail in ref. prb00, . The “quantum temperature” $`T_q`$ was shown to be given byprl2 ,

$$t=T_q/T_F=2/[1+0.86413(r_s^{1/6}1)^2]$$
(14)
where the Fermi temperature $`T_F`$ is simply the Fermi energy $`E_F=1/r_s^2`$ in our units. In effect $`T_q`$ is really a single-parameter representation of the density-functional 2D correlation energy.

The HNC equations for the PDFs of a classical fluid, and the Ornstein-Zernike(OZ) relations arehncref :

$`g_{ij}(r)`$ $`=`$ $`exp[\beta _{cf}\varphi _{ij}(r)+h_{ij}(r)c_{ij}(r)+B_{ij}(r)]`$
$`h_{ij}(r)`$ $`=`$ $`c_{ij}(r)+\mathrm{\Sigma }_sn_s{\displaystyle 𝑑𝐫^{}h_{i,s}(|𝐫𝐫^{}|)c_{s,j}(𝐫^{})}`$ (15)
This equation involves two quantities, viz., (i) the pair-potential $`\varphi _{ij}(r)`$, (ii) the bridge function $`B_{ij}(r)`$rosen ; yr2d . The other terms, e.g, $`c(r)`$, is the “direct correlation function” of the OZ equations, while $`h_{ij}(r)`$ = $`g_{ij}(r)`$-1.

Pair-potential. The $`\varphi _{ij}(r)`$ is the pair potential between the species $`i,j`$. For two electrons this is just the Coulomb potential $`V_{cou}(r)`$. If the spins are parallel, the Pauli principle prevents occupation of the same spatial orbital. This effect is equivalent to a “Pauli exclusion potential”, $`𝒫(r)`$ which when used in the HNC equation reproduces the $`g_{ij}^0(r)`$. Thus $`\varphi _{ij}(r)`$ becomes $`𝒫(r)\delta _{ij}+V_{cou}(r)`$. The Coulomb potential $`V_{cou}(r)`$ for a pair of point-charge electrons is $`1/r`$. However, an electron at the temperature $`T_{cf}`$ is localized to within a thermal wavelength. Thus, we use a “diffraction corrected” form

$$V_{cou}(r)=(1/r)[1e^{rk_{th}}]$$
(16)
In the case of the 2DEG $`k_{th}`$ was explicitely determined by numerically solving the Schrodinger equation for a pair of 2-D electrons in the potential $`1/r`$ and calculating the electron density in each normalized statepdwkth . That is, by solving the 2-D Schrodinger equation (only $`s`$-waves are needed for the $`r0`$ case):

$$\frac{d^2R}{dr^2}+\frac{1}{r}\frac{dR}{dr}+2m^{}(\frac{k^2}{2}\frac{1}{r})R=0$$
(17)
Here $`m^{}`$ = 1/2 is the effective mass of the electron pair. The radial function has the asymptotic form:

$$R(r,k)[\frac{2}{\pi kr}]^{1/2}cos(kr\frac{m^{}}{k}ln(r)+\delta _{2d})$$
(18)
The “on-top” density $`n(0)`$ for electrons at the temperature $`1/\beta `$ is:

$$\frac{n(0)}{n}=\frac{_0^{\mathrm{}}k^2𝑑kexp(\beta k^2/2m^{})R^2(0,k)}{_0^{\mathrm{}}k^d𝑑kexp(\beta k^2/2m^{})}$$
(19)
Using the same form for the diffraction correction, i.e., $`V_d(r)`$=exp$`(k_{th}r)/r`$, it was found that

$$k_{th}/k_{th}^0=1.158T_{cf}^{0.103}$$

where $`T_{cf}`$ is in a.u., and $`k_{th}^0`$ is the de Broglie thermal momentum $`(2\pi m^{}T_{cf})^{1/2}`$.

Bridge contribution. The “bridge” term arises from higher-order cluster interactions. They seem to play a role similar to the “back-flow corrections” used in QMC trial wavefunctionskwon . Since the Pauli exclusion already acts to reduce clustering of parallel spin electrons, the $`B_{ii}`$ bridge corrections are unimportant. In Ref. prl2  we included the $`B_{12}`$ terms via a simple hard-disk bridge function,yr2d  with a hard-disk packing fraction $`\eta _p=0.1175r_s^{1/3}(T_{cf}/E_F)`$. Since $`t=T_{cf}/E_F`$ is itself a function of $`r_s`$, the zero temperature packing fraction $`\eta _p(T=0)`$ is just a function of $`r_s`$. In generalizing the above expression to finite $`T`$, we note that $`r_s`$ is simply the coupling constant $`\gamma =P.E/K.E`$ per electron. We take the potential energy $`P.E`$ = $`1/r_s`$ for a p͡air of electrons, and evaluate the kinetic energy ($`K.E.`$) from the non-interacting internal energy as in Eq. 2. Then the finite temperature hard-disk packing fraction is given by

$$\eta _p=0.235\gamma ^{1/3}/[1+0.86413(\gamma ^{1/6}1)^2]$$
(20)
The next step is to use the $`\varphi _{ij}(r)`$, and solve the coupled HNC+Bridge and OZ equations for the two-component interacting  fluid. Evaluation at several values of the coupling constant and integration as in Eq. 13 leads to the correlation free energy.

Results. The fluid-phase with the lowest Helmholtz free energy $`F(r_s,T,\zeta )`$ is the stable phase. We have evaluated $`F`$ for the range $`t`$ =$`T/E_F<`$ 2, and 5 $`<r_s<40`$, and usually for five values of $`\zeta `$ and found that the stable phase is the unpolarized (paramagnetic) phase for higher temperatures and densities. The initial suggestion of Varsano and Senatore that the 2DEG at $`T`$=0 is spin polarized near $`r_s30`$, and the more recent discussion by Attaccalite et al., that a transition occurs at $`r_s`$=25.7 are greatly strengthened by our results at finite-$`T`$. At $`T`$=0 we find a $`\zeta =0\zeta =1`$ transition at $`r_s=r_c(\zeta :01)=26.7\pm 0.8`$, but the stabilization energy is within the inherent uncertainities of the CHNC method. However, when taken together with all the results for neighbouring finite-$`T`$ values, a more plausible picture emerges ( Fig. 1). When the temperature is increased, the value of $`r_c(\zeta :01)`$ decreases somewhat and then begins to increase. However, when $`t`$ reaches $`0.35`$, partially polarized states with $`\zeta 0.7`$ to 0.9 begin to compete for stability. This partially polarized state remains stable as long as the temperature is sufficiently low. An interesting aspect of the partially polarized regime is the existence of “ambi-spin” regions where the spin projection is NOT unambiguously determined since the energy differences as a function of the polarization are very small ( see below). For $`t>1.5`$, the stable phase becomes the usual unpolarized (paramagnetic) phase.

A plot of the energy difference $`\mathrm{\Delta }F(r_s,t,\zeta )=F(\zeta )F(0)`$, as a function of the spin polarization is shown in Fig. 2 at several temperatures and $`r_s`$ values close to one of the transition lines. Thus the $`t=0.05`$ curve is at $`r_s25.86`$ which corresponds to $`r_c(\zeta :01)`$. This curve is similar to that given in Fig. 1 of Ref. atta . However, the the transition regions near $`r_s26`$ and $`t0.35`$ provides an example of a phase where the Helmholtz fee energy is almost independent of the spin polarization. Such an “ambispin phase” would be very sensitive to any external effects (impurities etc.,) that we have not considered in this study. The regions in Fig. 1 close to the dotted line, begining from the ”triangular” region near $`r_s25`$ and $`t0.35`$ are such regions. The partially spin polarized state is quite unusual since, in some regions, a very light deviation from $`\zeta =1`$ leads to a significantly large stability. Thus the curve at $`t=0.5`$ is typical of densities near $`r_s3032`$ where the partially polarized stable is very close to unity but differs slightly from unity.

At $`T=0`$ the liquid phase becomes a Wigner crystal near $`r_s35`$. However, as the temperature is increased, this value will also increase. At low temperatures the crystal would be in equilibrium with the fully spin polarized phase. We have not investigated the equilibrium between the solid and liquid phases.

Thus the finite-$`T`$ 2DEG system provides many intriguing questions while complementing our understanding of the $`T=0`$ system. Our computer codes for generating the PDFs, the $`f_{xc}((r_s,T,\zeta )`$ and related properties of the 2D electron fluid may prove to be useful for further studies of this system. The codes may be accessed via the internetweb .

In conclusion, we have evaluated the total Helmholtz free energy $`F(r_s,T,\zeta )`$ = $`F_0+F_x+F_c`$, in the canonical ensemble, for a range of $`r_s<`$ 40, $`t2`$, and $`0z1`$. The terms $`F_0`$ and $`F_x`$ were evaluated by several methods, without resorting to perturbation expansions, while the $`F_c`$ was evaluated via a coupling constant integration of the pair-distribution functions obtained from the CHNC equations. The finite-$`T`$ results are consistent with recent claims, based on quantum Monte Carlo calculations, of the existence of a fully spin-polarized stable phase at $`T`$ =0 and $`r_s2526`$. Our results show that the finite-$`T`$ phase diagram of the clean, uniform 2DEG contains stable fully and partially spin-polarized phases. Ambispin phases where the free energy is quasi-independent of the spin polarization are also found. It is hoped that these results would stimulate new experiments and finite-$`T`$ QMC simulations.


text:
# Defect Structures in the Growth Kinetics of the Swift-Hohenberg Model

## I Introduction

What are the defects which control the long-time ordering of systems growing a striped pattern? This question arises in a variety of physical contexts. Here we are motivated by the recent experiments investigating the ordering of a two dimensional diblock copolymer system. The system studied offers a physical realization of the ordering in an isotropic two-dimensional smectic material. In these experiments they found that the late-time ordering satisfies scaling with a growth law $`Lt^x`$ with $`x=1/4`$ and the final stages of ordering are governed by the annihilation of sets of disclination quadrapoles. In this paper we address the question: Is the ordering in this physical system described by the Swift-Hohenberg (SH) model, the simplest model one can construct to govern the ordering in stripe forming systems?

We investigate the growth kinetics of the Swift-Hohenberg model for a small control parameter ($`ϵ=0.1`$) in two dimensions and quenches to zero temperature. It is this regime which appears most likely to correspond to the experimental situation. In large $`ϵ`$ regime the system evolves to a glassy state. We focus primarily on the defect structures generated in the ordering of the system. In the most naive picture of this ordering process one can think in terms of an initial local layering, as in a smectic, in some direction. This ordering can be disrupted by point defects: dislocations and disclinations. This suggests a coarsening picture with annihilating point defects similar to the case of the XY model and a growth law with exponent $`x=1/2`$. This simple picture is not seen in simulations. We find, in agreement with the numerical results of Hou et al.  and Boyer and Viñals , that the defect structures for the SH model are dominated by grain boundaries which persist for long times. Unlike the case of an XY model, the ordering is not dominated by annihilation of isolated point defects. These are observed but are not the dominant structures.

We find numerically, at late times after finite size effects enter, that the system becomes anisotropic, and the grain boundaries shrink. In this case one sees a cross-over to an effective growth exponent $`x=1/2`$.

We give below a detailed numerical study of the statistical properties of the defects disrupting striped pattern formation in the SH model. In order to carefully discuss the defects we need a reliable filter for finding them. We present an algorithm which effectively locates defects and grain boundaries for any control parameter $`ϵ`$. We can distinguish between grain boundaries and other defects, and track their trajectories separately. We compare this method to the other approaches used in earlier work in appendix A.

There are a number of ways of characterizing the degree of ordering in these systems: (i) Counting the number and size of defects and their evolution with time. (ii) Monitoring the lowering of the average effective driving free energy as a function of time. (iii) Evaluation of the nematic order parameter correlation function and its associated scaling behavior. (iv) Evaluation of the order parameter structure factor and its associated scaling behavior. We find that (i), (ii), and (iii) can all be characterized by a single growth law with the exponent near $`1/3`$, while the order parameter scaling, as found by others, is characterized by a growth law with the exponent near $`1/5`$.

## II Swift-Hohenberg Model

The Swift-Hohenberg model for a scalar order parameter, $`\psi `$, is specified by the equation of motion

$`{\displaystyle \frac{\psi }{t}}=ϵ\psi \psi ^3\left(q_0^2+^2\right)^2\psi +\zeta ,`$ (1)
where $`ϵ`$ is a positive control parameter, $`q_0`$ is the magnitude of an ordering wavenumber and $`\zeta `$ is the Gaussian noise satisfying $`\zeta (𝐫,t)\zeta (𝐫^{},t^{})=2\mathrm{\Gamma }\delta (𝐫𝐫^{})\delta (tt^{})`$, where the noise strength $`\mathrm{\Gamma }`$ is proportional to the final temperature governing the system after a quench. We will focus here on quenches to zero temperature where we can set $`\mathrm{\Gamma }`$ and the noise $`\zeta `$ to zero. We are interested in the growth kinetics problem where we prepare this system initially in a completely disordered state. We then allow the system to evolve forward in time to form a striped pattern. For example one could choose

$`\psi (𝐱,t_0)\psi (𝐲,t_0)=\mathrm{\Psi }_0^2\delta \left(𝐱𝐲\right),`$ (2)
where $`\mathrm{\Psi }_0^2`$ is a constant. However the precise form of the initial conditions is not important .

This model can be formulated as a Langevin equation driven by an effective Hamiltonian:

$`_E={\displaystyle d^2x\left\{\frac{ϵ}{2}\psi ^2+\frac{1}{4}\psi ^4+\frac{1}{2}\left[\left(q_0^2+^2\right)\psi \right]^2\right\}}.`$ (3)
If we introduce

$`E(t)=_E_t,`$ (4)
where the average is over an ensemble of initial conditions, then $`E(t)`$ is lowered as the system orders in a striped pattern with wavenumber $`q_0`$.

Eventually the system approaches an ordered state described approximately by the single-mode approximation  where, assuming layering along the z-direction,

$`\psi _0=A_0\mathrm{cos}q_0z.`$ (5)
If we put this ansatz into Eq.(3), assume that the system is an integral number of wavelengths in the z-direction, and minimize with respect to the amplitude $`A_0`$, we obtain the results,

$`A_0^2`$ $`=`$ $`{\displaystyle \frac{4ϵ}{3}},`$ (6)
$`\psi _0^2`$ $`=`$ $`{\displaystyle \frac{q_0}{2\pi }}{\displaystyle _0^{2\pi /q_0}}(A_0\mathrm{cos}q_0z)^2𝑑z={\displaystyle \frac{2ϵ}{3}},`$ (7)
and

$`E_{eq}={\displaystyle \frac{ϵ^2}{6}}S,`$ (8)
where $`S`$ is the area of the system. Pomeau and Manneville have shown that this is a very good approximation for the “ground” state even for moderately large values of $`ϵ`$. In the growth kinetics context the approach to equilibrium is monitored by

$`\mathrm{\Delta }E(t)E(t)E_{eq}L_E^1(t),`$ (9)
and

$`\mathrm{\Delta }\psi ^2(t)\psi _0^2\psi ^2_tL_\psi ^1(t),`$ (10)
where $`L_E(t)L_\psi (t)`$ .

Another measure of the ordering in the system is given by considering the director field

$`\widehat{n}(𝐱)={\displaystyle \frac{\psi (𝐱)}{|\psi (𝐱)|}},`$ (11)
and the associated nematic order parameter

$`Q_{\alpha \beta }=Q_0\left[\widehat{n}_\alpha \widehat{n}_\beta {\displaystyle \frac{1}{2}}\delta _{\alpha \beta }\right].`$ (12)
In two dimensions, however, all of the information in this order parameter is contained in the quantity $`\mathrm{cos}2\theta `$ where $`\widehat{n}=(\mathrm{cos}\theta ,\mathrm{sin}\theta )`$. It is easy to show, for example, that

$`C_{nn}(𝐱,𝐲,t)2\mathrm{Tr}Q(𝐱,t)Q(𝐲,t)_t`$ (13)
$`=\mathrm{cos}[(\phi (𝐱,t)\phi (𝐲,t)]_t.`$ (14)
where

$`\phi (𝐱,t)=2\theta (𝐱,t).`$ (15)
If we define

$`\widehat{B}_x=\widehat{n}_x^2\widehat{n}_y^2`$ (16)
$`\widehat{B}_y=2\widehat{n}_x\widehat{n}_y`$ (17)
then

$`C_{nn}(𝐱,𝐲,t)=\widehat{B}(𝐱,t)\widehat{B}(𝐲,t)_t.`$ (18)
The nematic order parameter correlation function, $`C_{nn}`$, was shown by Christensen and Bray  to obey scaling in the conventional form

$`C_{nn}(𝐫,t)=F\left(r/L_n(t)\right),`$ (19)
where $`𝐫=𝐱𝐲`$. Elder, Viñals and Grant  showed that the scaling of the order parameter structure factor

$`S(k,t)=\left|\psi _𝐤(t)\right|^2=L_s(t)F_1((kq_0)L_s(t)),`$ (20)
differs from that observed in ordering system without stripes: $`S(k,t)=L^2(t)F_2(kL(t))`$.

## III Review of Previous Work

The early work on this problem focused on establishing the final equilibrium state reached after a quench. This is a two dimensional system and by forming stripes one has a broken continuous symmetry. The behavior of the system at non-zero temperatures, as for the two dimensional X-Y model, requires, as pointed out by Toner and Nelson  a treatment of both long wavelength fluctuations in the layers and free defects. Above a Kosterlitz-Thouless type transition one has an isotropic phase while below this transition one has a phase with persistent orientational order.

In an early paper, Elder, Viñals, and Grant carried out a numerical analysis leading to the scaling solution given by Eq.(19). Working with fixed $`ϵ=0.25`$ they looked at the system’s ordering as a function of noise strength $`\mathrm{\Gamma }`$. They found a qualitative difference between low noise and high noise. For the large noise case they found a rapid (exponential) relaxation to the asymptotic stationary state and a power-law approach for the lower noise case. Their results are in agreement with the picture due to Toner and Nelson that one has a transition to an isotropic state for large enough noise. There is no real ordering in the isotropic state and this is why there is exponential decay to the equilibrium state. In the ordered state one has scaling and a power-law growth law which, for small noise, they found to have an exponent $`x_s=1/4`$. They found a smaller exponent $`x_s=1/5`$ at low temperatures, but they had less statistics and there appeared to be ”difficulty removing defects”. They argued for a late time cross over to the expected $`x=1/2`$ but they did not see this.

Cross and Meiron  also studied the SH model numerically in the absence of noise. They found a $`x_s=1/4`$ for $`ϵ=0.25`$. The dynamics appear to freeze for higher $`ϵ`$. They looked at the defect structure but in a qualitative way noting the existence of domain walls rather than a set of isolated point defects. The theoretical discussion in their paper is based on the phase-field approximation

$`{\displaystyle \frac{\varphi }{t}}=\left(D_{}_{}^2+D_{}_{}^2\right)\varphi ,`$ (21)
which from the most naive point of view suggests a growth law with exponent $`x=1/2`$. They discuss some selection mechanisms which could lead $`D_{}`$ and $`D_{}`$ to adjust themselves to zero and reduce $`x`$ to $`1/3`$ or $`1/4`$. They concluded that they did “not have a good theoretical understanding of these results” and suggested that the defects in the problems should be treated explicitly.

Hou, Sasa, and Goldenfeld (HSG)  confirmed previous numerical results which showed for $`ϵ=0.25`$, $`x_s=1/5`$ with zero noise and $`x_s=1/4`$ with nonzero noise as obtained from the structure factor scaling. They went further and used a simple method to identify domain walls and measure their lengths (more about this below). They measured excess energy, $`\mathrm{\Delta }E(t)`$, and the domain wall length and found that they show the same scaling exponents $`1/4`$ at zero noise and $`0.3`$ at non-zero noise. The energy does go to the lowest order in $`ϵ`$ value of $`ϵ^2/6`$ in the noiseless limit. They find “defects are indeed the driving force behind the coarsening process due to its dominant contribution to the excess energy.” They suggest that the phase field approach gives the wrong exponent because it does not include the effects of defects. For larger $`ϵ`$ (=0.75) they found much slower logarithmic growth. The system seems to become glassy.

Christensen and Bray  also carried out numerical work on the SH model for $`ϵ=0.25`$ and found $`x_s=1/5`$ for zero noise and $`x_s=1/4`$ for nonzero noise. From scaling of the director correlation function they find exponents are $`0.25`$ and $`0.30`$ for zero and nonzero noise. They suggest that there is a cross over to $`x=1/2`$ at very long times. The theory they developed does not include defects.

Boyer and Viñals  point out ”Near the bifurcation threshold, the evolution of disordered configurations is dominated by grain boundary motion through a background of largely immobile curved stripes”. They find for small $`ϵ`$ an exponent $`x=1/3`$ which they interpret as arising from a law of grain boundary motion . Elsewhere  they also point out for larger values of $`ϵ`$ the dynamics cross over to a frozen state with quenches to zero temperature. This glassy behavior is associated with grain boundary pinning.

## IV Numerical Results for SH Model

We present here our numerical results for the SH equation. We follow the numerical prescriptions of Bray and Christensen . We use the finite difference scheme on two dimensional lattice of sizes $`256\times 256`$ and $`512\times 512`$ with periodic boundary conditions. We set $`ϵ=0.1`$, $`\mathrm{\Delta }r=\pi /4`$ and $`\mathrm{\Delta }t=0.03`$. We replace $`_t\psi (𝐫,t)`$ by $`\left(\psi _{ij}^{n+1}\psi _{ij}^n\right)/\mathrm{\Delta }t`$, and $`^2\psi (𝐫,t)`$ by

$$^2\psi _{ij}=\frac{1}{(\mathrm{\Delta }x)^2}\left[\frac{2}{3}\underset{NN}{}+\frac{1}{6}\underset{NNN}{}\frac{10}{3}\right]\psi _{ij}$$
(22)
where $`NN`$ and $`NNN`$ mean the nearest neighbors and next-nearest neighbors respectively. By choosing the proper scale of time, space and the order parameter, we can set $`q_0=1`$. The systems have eight grid points per wavelength. We used uniformly distributed random initial conditions.

For the smaller $`256\times 256`$ systems we were able to follow the ordering process to very late stages. Some of the independent trials proceed to a final state where we have a set of well aligned layers.

In Fig. 1 we plot $`\mathrm{\Delta }E(t)`$ and $`\mathrm{\Delta }\psi ^2(t)`$ for an ensemble of runs on a $`256\times 256`$ lattice. We note that there are two regimes where $`L_E`$ defined by Eq.(9) is described by different exponents. For $`t_s<t<t_c`$ ($`t_s300`$ and $`t_c9000`$) we find $`x_E0.3`$, while for $`t>t_c`$ we find $`x_E0.5`$. The cross over at $`tt_c`$ appears to be due to the finite size effects, as we discuss below. For $`t>t_c`$ the system is effectively anisotropic and we find an effective exponent $`x_E`$ near to $`1/2`$.

In Fig. 2 we plot $`\mathrm{\Delta }E(t)`$ and $`\mathrm{\Delta }\psi ^2(t)`$ for a $`512\times 512`$ system. In this case we see that $`t_c`$ has been extended to much larger values and we have not been able to follow the ordering process to completion. Our fits to $`\mathrm{\Delta }E(t)`$ and $`\mathrm{\Delta }\psi ^2(t)`$ in the regime $`t_s<t<t_c`$ again gives, to higher accuracy, $`x_E=x_\psi =1/3`$.

To probe directly the stripes’ increasingly orientational order, we measure the nematic order parameter correlation function $`C_{nn}(𝐫,t)`$ in the $`512\times 512`$ system. The results, averaged over $`57`$ runs, are shown in Fig. 3. We obtain scaling with a correlation length obeying the growth law $`L_nt^{0.36}`$. We can estimate the time $`t_c`$ when the cross-over begins in this larger system as follows. The system becomes anisotropic and one expects cross-over when the correlation length $`L_n`$ grows to be some substantial fraction of a lateral dimension of the system. In terms of ratios we can write

$$\frac{L_n(t_c(512))}{L_n(t_c(256))}\frac{512}{256}=\left[\frac{t_c(512)}{t_c(256)}\right]^{1/3}.$$
(23)
In the $`256\times 256`$ system $`t_c(256)9000`$, so we obtain $`t_c(512)60000`$. Notice that in Fig. 2 the effective exponent $`x_E`$ begins to increase at the time $`5000070000`$, which is consistent with our estimate.

In Fig. 4 we plot the structure factor $`S(𝐤,t)`$ and show that scaling holds in the form given by Eq.(20) with a growth law characterized by an exponent $`x_s=0.24`$ as shown in the insert. Our results here agree with those found previously that the exponent governing the growth law for the structure factor is significantly smaller than that governing the nematic order parameter.

## V Defect Structures and dynamics

In Fig. 5 we show a typical configuration for the Swift-Hohenberg model for a quench to zero temperature after a time $`12000`$ for a $`512\times 512`$ system. Notice the rather complicated structure which includes dislocations, disclinations and grain (domain) boundaries. Our main focus in this paper is to study the statistics of these defects. In appendix A we discuss an algorithm for picking out the defects and tracking their motions.

If we look at Fig. 5 we see that it shows a complicated situation with a variety of different defect structures which one can identify by eye at the length scale of several layer spacings. At a more fundamental level we need a way of identifying which points in space, at the level of each site on the numerical grid, are part of a defect. At the shortest length scale in the problem the order parameter is $`Q_{\alpha \beta }`$ defined by Eq.(12). For this two-dimensional system this can be replaced by the vector order parameter $`\widehat{𝐁}`$ defined by Eqs.(16) and (17). The assumption is that all of the defects in the system can be built up from the $`\pm \frac{1}{2}`$ disclinations in the director field $`\widehat{n}`$ which translate into vortices with charge $`\pm 1`$ for the field $`\widehat{𝐁}`$. We identify these defects by looking for the cores of the vortices. We can find the cores of the defects by looking from those sites where $`\widehat{𝐁}`$ is changing rapidly. We can define

$`A={\displaystyle \underset{\alpha ,\beta }{}}\left(_\alpha B_\beta \right)^2`$ (24)
and identify defect points as those sites where $`A`$ is larger than some value. Notice that $`A`$ can also be written in the form

$`A=4{\displaystyle \underset{\alpha ,\beta }{}}\left(_\alpha n_\beta \right)^2=\left(_\alpha \phi \right)^2.`$ (25)
The precise numerical determination of $`A`$ is discussed in appendix A. Notice that $`A`$ is proportional to the gradient energy for an isotropic nematic.

In analyzing their experimental data Harrison, et al.  found a set of fundamental disclinations and from these built up dislocations as bound disclinations with opposite charge. They used this procedure to identify a large dislocation density. Most of the fundamental disclinations went into forming these dislocations since in the end the ratio of dislocations to the remaining disclinations was about ten to one. In our case the situation is complicated by the grain boundaries. We first separate the defects into compact point defects and larger grain boundaries. For the point defects we determine the topological charge by taking the usual phase-angle path integral around the center of mass of the defect. Those defects with plus or minus unit charge are identified as disclinations, while those with zero charge are dislocations. Then we can track the motion of each single defect.

As an example of the method we show in Fig. 6 the path of a dislocation. We see that some dislocations travel over long distances during very long times. It seems that the dislocations are more stable when compared to grain boundaries and disclinations. Our simulations on $`256\times 256`$ system show that after the annihilation of point defects and grain boundaries, some dislocations still exist in the system. Our simulations show that there are also dislocations which are pinned and move little.

The number of disclinatins is quite small. And we notice that they are rather immobile, which is consistent with Boyer and Viñals’ discussion .

The most important motion of grain boundaries is that they can move over long distances and combine with other grain boundaries. As shown in Fig. 7, two grain boundaries can combine to form a larger grain boundary. Thus the number of grain boundaries decreases while their average size increases. This process happens on a time scale of the order $`1000`$ dimensionless time units.

As shown in Fig. 8, one grain boundary can sweep across a quite large area. At the same time its size decreases. This process occurs on a time scale of the order $`10000`$. According to our observations, the grain boundaries’ motions also relieve the stripe curvatures through disclination annihilations. After one grain boundary passes through a disclination, the disclination disappears. This is consistent with Boyer and Viñals’ prediction .

Next we focus on the statistics of the defects generated by the model. In Fig. 9 we plot the total number of points in grain boundaries, dislocations and disclinations separately for the $`256\times 256`$ system. We see that the grain boundaries dominate. In the scaling regime ($`t_s<t<t_c`$), we see that the number of points corresponding to grain boundaries and all defect points, the curves $`a`$ and $`b`$ can be fit to $`t^{1/3}`$. At late stages the disclinations disappear, while the dislocations and grain boundaries persist. The number of disclinations decreases much faster than the other defects.

In Fig. 10 we plot the average number of grain boundaries $`\overline{n}`$ and the average size of a grain boundary $`\overline{l}`$ for the $`256\times 256`$ system. We use the number of points in one grain boundary as a measure of its size. For $`t_s<t<t_c`$, $`\overline{n}`$ decreases but $`\overline{l}`$ increases due to the combining of grain boundaries. The shrinkage of their sizes is not as important as the combinations. However for $`t_c9000`$ the correlation length $`L_n`$ is the same order as the system’s size, and the large grain boundaries stop growing. After that the shrinkage is important . In the scaling regime ($`t_s<t<t_c`$), $`\overline{n}t^{0.45}`$, and $`\overline{l}t^{0.13}`$. So $`\overline{n}\overline{l}t^{1/3}`$, which is consistent with Fig. 9.

In Fig. 11 we plot the total number of points in grain boundaries, dislocations and disclinations separately for the $`512\times 512`$ system. The number of points in disclinations is much smaller than that of dislocations and at very late stages it decreases to the order of 1, which in fact indicates the disappearance of disclinations. Now the scaling regime extends to much longer times. The domain walls and dislocations’ scaling exponents are both $`1/3`$, which is same to the scaling of the energy. However, the scaling exponent of disclinations, given by 0.57 is much larger. So we conclude that disclinations are not the dominant structures in SH system.

In Fig. 12 we plot the number of grain boundaries $`\overline{n}`$ and the average size of a grain boundary $`\overline{l}`$ for $`512\times 512`$ system. All of our data falls in the scaling regime.. The plot of the average number $`\overline{n}`$ of grain boundaries versus time $`t`$, can be fit to $`\overline{n}t^{0.49}`$ and the average size $`\overline{l}`$ v.s. time $`t`$, can be fit by $`\overline{l}t^{0.17}`$. So we have $`\overline{n}\overline{l}L^1t^{1/3}`$, which is consistent with the results shown in Fig. 11.

Although we did not count the number of disclinations directly, it is proportional to the number of lattice points in disclinations, i.e. $`t^{0.57}`$. This is because the average number of lattice points in one point disclination, which is about $`1020`$, is quite stable during the simulation. By the same reasoning, we find that the number of dislocations is proportional to $`t^{1/3}`$. It is interesting to note that the number of grain boundaries scales as $`t^{0.49}`$. This exponent is near to that for disclinations.

The number of grain boundaries is about 5 at $`t70000`$, the number of disclinations is about $`0`$ or $`1`$ at $`t50000`$, and the number of dislocations is on the order of 10 at $`t50000`$, as can be seen in Figs. 11 and 12.

Since we can track the motion of each defect, we can measure their speeds. We define the speed of each as the speed of its mass center (see appendix A for the definition of “mass center”). If in a time $`\mathrm{\Delta }\tau `$, the mass center travels over a distance $`\mathrm{\Delta }d`$, then the speed is $`v=\mathrm{\Delta }d/\mathrm{\Delta }\tau `$. If $`\mathrm{\Delta }\tau `$ is small enough, we found that $`v=0`$ has the biggest probability.. If $`\mathrm{\Delta }\tau `$ is large enough, all the details are coarse-grained and we observe a continuous distribution of the speed and the largest probability appears at a non-zero speed, as is shown in Fig. 13 where $`\mathrm{\Delta }\tau =60`$. We measured the speed distributions of domain walls and point defects separately. As we have already seen, for point defects the number of disclinations is much smaller than that of dislocations, so what we measure in the latter case is in fact the speed distribution of dislocations.

The speed distribution has a long tail which decreases as a power law. The numerical fits at different times give us different exponents. However the tail exponents at different times do distribute in a narrow region, as is shown in Fig. 14. The exponents of grain boundaries are quite different from those of point defects. If we ignore the exponents at very early times when grain boundaries just begin to form and at the very late times when the grain boundaries have already disappeared, the mean value of the grain boundaries’ exponents is $`1.50`$, and that of point defects’ exponents is $`2.10`$.

We also measured the average speed of the point defects and grain boundaries as a function of time after the quench, as is shown in Fig. 15. The average speed of point defects decreases as $`t^{0.48}`$; the average speed of grain boundaries goes as $`t^{0.35}`$. The scattering of the points at late stages is due to the small data sample at those times..

## VI Summary

We have studied the dynamics of the defect structures in the SH model after quenches to zero temperature with a control parameter of $`ϵ=0.1`$. We find in agreement with earlier workers that the kinetics in the ordering regime, before finite-size effects enter, are dominated by the existence of moving and coalescing grain boundaries. In this regime the average size of these grain boundaries is growing and they are relatively mobile. Under the influence of finite size effects these grain boundaries shrink, the system becomes anisotropic and the ordering process speeds up.

We also measured the speed distribution of all structures that appear in the system. The average speed is decreasing as a power law and the distributions show a power-law behaviors at large speeds. However we can only get a rough estimate due to the poor statistics.

Let us return to the question of whether the SH model gives a good description of the physical system studied by Harrison, et al.. The SH model, for small control parameter $`ϵ`$, does give coarsening with an exponent in roughly the same range as in the experiment ($`1/41/3`$). The ordering is constrained to be slower then the picture where one has a simple point defect pair annihilation process. However, the defect structures in the SH model and experiment appear quite different. The disclination quadrapole annihilations seen in experiments are not observed in the late stages of the evolution of the SH model. In the SH model grain boundaries dominate the evolution in the scaling regime, but these structures appear to play a limited role in the experiments. We must make clear that our numerical results are for systems with many fewer roll periods compared to the experimental systems ($`10^2`$ compared to $`10^5`$), so it is possible that things change as we increase the size of the ordering system. However, our study of 256 and 512 systems shows that they differ only in the time when the finite-size effect enters. This indicates that a even larger system will display the same behavior except that the finite-size effect enters at a even later time. So we conclude that the SH model does not give a physically faithful description of the ordering in the experimental system. This raises the provocative question: Are there many different types of scenarios for ordering striped systems? We will address this question by looking at other competing models for striped formation elsewhere.

Acknowledgments: We thank Dr. C. Harrison for providing us with a copy of his thesis. This work was supported by the National Science Foundation under Contract NO. DMR-0099324.

## A The algorithm for picking out defects and finding domain walls

Hou, et al.  proposed a method, the HSG method, to measure the length of grain boundaries. They computed the quantity $`A^2\psi ^2+(\psi )^2/q_0^2`$, and if the calculated $`A^2`$ is bigger than an upper filter $`0.7\times Max(A^2)+0.3\times AVG(A^2)`$ or smaller than a lower filter $`0.7\times Min(A^2)+0.3\times AVG(A^2)`$, that point is counted as belonging to a domain wall. When $`ϵ`$ is small, this method gives quite good results.

However, when $`ϵ`$ increases, the original filters are no longer applicable. They fail to pick out most of the points and the filters must be re-chosen. For example, at $`ϵ=0.6`$, the filters $`0.5\times Max(A^2)+0.5\times AVG(A^2)`$ and $`0.5\times Min(A^2)+0.5\times AVG(A^2)`$ can give a satisfying result; while for $`ϵ=0.75`$, $`0.4\times Max(A^2)+0.6\times AVG(A^2)`$ and $`0.4\times Min(A^2)+0.6\times AVG(A^2)`$ are the better choices. Sometimes this method is unable to pick out all the defects for any choice of filter.

We introduce here a method which works for all $`ϵ`$ and picks out all of the defects and nothing more.

First let us define some useful quantities. Suppose the system is discrete on the $`xy`$ plane, with $`𝐱=(i,j)`$ (square lattice). At a fixed time, starting from the order parameter field $`\psi (𝐱)=\psi _{i,j}`$, we can define a director field $`\widehat{n}(𝐱)`$ as given by Eq.(11) where $`\psi (𝐱)`$ is defined by the usual finite difference scheme, i.e.

$`\psi (𝐱)=({\displaystyle \frac{\psi _{i+1,j}\psi _{i1,j}}{2\mathrm{\Delta }r}},{\displaystyle \frac{\psi _{i,j+1}\psi _{i,j1}}{2\mathrm{\Delta }r}}),`$ (A1)
where $`\mathrm{\Delta }r`$ is the lattice space of the system. In two-dimensional cases the nematic order parameter, $`Q_{\alpha \beta }`$ is completely specified by the angle

$$\phi (𝐱)=2\theta (𝐱),$$
(A2)
where

$$\theta (𝐱)=\mathrm{arctan}\left(\frac{\widehat{n}_y(𝐱)}{\widehat{n}_x(𝐱)}\right).$$
(A3)
Rather than using $`\phi (𝐱)`$ given by the two equations above we introduce some local smoothing. First we compute

$`\widehat{B}_y=\mathrm{sin}\phi (𝐱)=2\widehat{n}_x(𝐱)\widehat{n}_y(𝐱),`$ (A4)
$`\widehat{B}_x=\mathrm{cos}\phi (𝐱)=2\widehat{n}_x(𝐱)^21.`$ (A5)
Then we smooth these two fields using the iterative process:

$$f_{(n+1)}(i,j)=\frac{1}{2}f_{(n)}(i,j)+\frac{1}{8}\underset{(i^{},j^{})NN}{}f_{(n)}(i^{},j^{}),$$
(A6)
where $`f_{(n)}`$ is $`\mathrm{sin}\phi `$ or $`\mathrm{cos}\phi `$ after $`n`$ iterations, and $`NN`$ means the $`4`$ nearest neighbors of $`(i,j)`$ on the square lattice. This process will suppress the small fluctuations of $`\phi (𝐱)`$ away from the defects, while the variation of of $`\phi (𝐱)`$ near a defect core remains large. Our calculations show that $`5`$ iterations provides a sufficiently smooth set of fields for our purposes. In the next step, we calculate $`\phi (𝐱)`$ from $`\mathrm{sin}\phi (𝐱)`$ and $`\mathrm{cos}\phi (𝐱)`$ using

$$\phi (𝐱)=\mathrm{arctan}\left[\frac{\mathrm{sin}\phi (𝐱)}{\mathrm{cos}\phi (𝐱)}\right],$$
(A7)
where we adopt the convention that $`\pi <\phi (𝐱)<\pi `$.

In picking a filter we want to look at the spatial variation of the $`\phi (𝐱)`$ field. $`\phi (𝐱)`$ can be evaluated as for $`\psi (𝐱)`$. However, there is a subtlety here. For example, if $`\phi _{i+1,j}=\pi \delta \varphi _1`$ and $`\phi _{i1,j}=\pi +\delta \varphi _2`$, where $`\delta \varphi _1`$ and $`\delta \varphi _2`$ are small angles, then the difference between the nematic tensor $`Q_{\alpha \beta }(i+1,j)`$ and $`Q_{\alpha \beta }(i1,j)`$ should be a small quantity. But $`(\phi _{i+1,j}\phi _{i1,j})/2=\pi (\delta \varphi _1+\delta \varphi _2)/2\pi `$, which means that if we calculate the change rate of $`\phi (𝐱)`$ in exactly the same way as Eq.(A1), we will get a wrong answer in this context. To avoid such a problem, we define the difference between $`\phi (𝐱)`$ and $`\phi (𝐱^{})`$ as the quantity with the smallest absolute value among the choices $`\phi (𝐱^{})\phi (𝐱)`$ and $`\phi (𝐱^{})\phi (𝐱)\pm 2\pi `$. And we use this quantity in determining

$`A(𝐱)=|\phi (𝐱)|^2`$ (A8)
which is the key quantity in our analysis.

Our method is based on the observation that at the core region of a defect (dislocation, disclination or part of a domain wall) the angle field $`\phi (𝐱)`$ changes rapidly, while in the region away from the defect’s core the $`\phi (𝐱)`$ field is rather smooth, as can be seen in Fig. 16. Thus we can conclude with confidence that those points with larger change rates of $`\phi (𝐱)`$ must belong to some defect’s core region or a part of a grain boundary.

We find that $`A(𝐱)0`$ away from defects, but increases very rapidly in the vicinity of any defect. An example is given in Fig. 17. Therefore as long as $`A(𝐱)`$ is large enough, we can identify the point $`𝐱=(i,j)`$ as part of the core region of a defect. Naturally we set up a threshold $`A_0`$, and any point with $`4(\mathrm{\Delta }r)^2A(𝐱)>A_0`$ is counted as belonging to some defect’s core. Because the value of $`A(𝐱)`$ is much larger in the defects’ cores than at any other places, a range of values of $`A_0`$ can be used to find the positions of the defects core regions. With a smaller threshold the program will pick out more points in the core regions, and with a larger one it will pick out fewer points in the core regions. Our experience shows that if $`A_0`$ takes the value of $`210`$, the program picks out the same defects cores and grain boundaries. As is shown in Fig. 18 and Fig. 19, it picks out all the defects without irrelevant points.

After we have used the above algorithm to pick out the points in the core regions of the point defects and grain boundaries, we can distinguish between these two structures. The difference between them is obvious. The point defects are compact in space while the grain boundaries are ramified.

First, we must group the points we have identified according to whether or not they are in the same structure. The points in one point defect core or grain boundary are picked out because the director field changes drastically on those sites. They are very near to each other. However, they may not be neighbors. So we define a filter $`a_0`$, and when any two points’ distance is less than $`a_0`$, they are supposed to be in the same defect or grain boundary’s core region. We use the cluster multiple labeling method of Hoshen and Kopelman  to pick out such point clusters. Thus given the system’s status at any time, we can find those sets of points corresponding to each individual defect or grain boundary.

Now we measure the approximate size of these structures and then distinguish between point defects and grain boundaries. We use the number of the points in the set as the size of the corresponding structure. This approximation reflects the actual size of the defect or grain boundary quite well. Then we define a filter $`l_0`$, and when the structure’s size is larger than $`l_0`$, we regard the corresponding structure as a grain boundary, otherwise it’s taken to be a point defect (dislocation or disclination).

We employed $`a_0=5\mathrm{\Delta }r`$ where $`\mathrm{\Delta }r`$ is the lattice spacing and $`l_0=18`$. The results are quite satisfying.

After we have picked out the point defects, we can devide them into disclinations and dislocations. We follow Harrison’s method . Given the angle field $`\phi (𝐱)`$ computed in Eq.(A6) after the smoothing process, we do an integral of the variation of $`\phi (𝐱)`$ over a counterclockwise close path around the “mass-center” of a point defect, which is defined below. The condition for a defect to be a disclination is

$$\frac{\phi }{s}𝑑s=\pm 2\pi $$
(A9)
The integral is zero if the defect is a dislocation. To make the computation easier, we choose a $`16\times 16`$ square with the mass-center at its center as the integration route.

To record the motion of one single defect or grain boundary, we track the motion of the corresponding point set’s “mass center”, which is defined as follows. Suppose the point set has $`n`$ points with coordinates $`𝐫_i,i=1,2,\mathrm{},n`$. Then the “mass center” of the point set is defined as $`𝐫={\displaystyle \underset{i}{}}𝐫_i/n`$, just like the usual mass center definition in classical mechanics but with all masses equal to one.

In the evolution of the SH model, we sample the system every $`500`$ time steps (in our case this is equal to 15 dimensionless time units), which is a quite short-time period in the simulation. We then identify all the dislocations, disclinations and grain boundaries, distinguish among them, and compute their centers of mass. Suppose at time $`t_1`$, we have the set of mass centers $`P=\{𝐩_i,i=1,2,\mathrm{},n_1\}`$, and at time $`t_2`$, the “mass center” set is $`Q=\{𝐪_j,j=1,2,\mathrm{},n_2\}`$; usually $`n_1n_2`$. Define $`d_{PQ}(i,j)=|𝐩_i𝐪_j|`$. We assume that the defects and grain boundaries do not move much in such a short time period. So if there exist two integers $`k[1,n_1]`$ and $`l[1,n_2]`$, such that

$$d_{PQ}(k,l)=\underset{j[1,n_2]}{\mathrm{min}}d_{PQ}(k,j)=\underset{i[1,n_1]}{\mathrm{min}}d_{PQ}(i,l),$$
(A10)
it is quite reasonable to believe that $`𝐩_k`$ and $`𝐪_l`$ are just the same defect’s or grain boundary’s “mass center” at two successive times. Using this method, we are able to find out the trajectories of the “mass centers” as time goes on. Not all points in $`P`$ and $`Q`$ can be grouped into such pairs. On the one hand, this is because $`n_1n_2`$; on the other hand, this is also due to the criterion (A9) applied onto $`𝐩_k`$ and $`𝐪_l`$. Physically, this is consistent with the phenomena of the defect annihilation and the combination, split and shrinkage of grain boundaries.

## B Measurement of the Nematic Correlation Function

To probe the stripes’ increasingly orientational order, we define the correlation function which is similar to the one employed by Christensen and Bray.

$`C_{nn}(𝐫,t)={\displaystyle \frac{1}{N^2}}{\displaystyle \underset{𝐱}{}}\mathrm{cos}[\phi (𝐱+𝐫,t)\phi (𝐱,t)]`$ (B3)
$`={\displaystyle \frac{1}{N^2}}{\displaystyle \underset{𝐱}{}}\mathrm{cos}\phi (𝐱+𝐫,t)\mathrm{cos}\phi (𝐱,t)+`$
$`+{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{𝐱}{}}\mathrm{sin}\phi (𝐱+𝐫,t)\mathrm{sin}(𝐱,t),`$
where $`N^2`$ is the area of the system and the angular brackets denote the statistical average over different initial conditions. The definition of the angle $`\phi (𝐱)`$ is given in appendix A.

Now in Eq.(B3) the function has been split into two parts which have the same form

$$G(𝐫,t)=\frac{1}{N^2}\underset{𝐱}{}f(𝐱+𝐫,t)f(𝐱,t),$$
(B4)
with $`f(𝐱,t)=\mathrm{cos}\phi (𝐱,t)`$ and $`f(𝐱,t)=\mathrm{sin}\phi (𝐱,t)`$ separately. Eq.(B4) can be easily calculated by fast Fourier transformation (FFT). First FFT $`f(𝐫,t)`$ to obtain its Fourier components $`\stackrel{~}{f}(𝐤,t)`$. Then $`\stackrel{~}{G}(𝐤,t)=|\stackrel{~}{f}(𝐤,t)|^2`$, and inverse Fourier transformation gives $`G(𝐫,t)`$. We compute the two parts in Eq.(B4) separately, and then add to obtain the correlation function $`C_{nn}(𝐫,t)`$.