text:
# Magnetic-field dependence of electron spin relaxation in n-type semiconductors

## Abstract

We present a theoretical investigation of the magnetic field dependence of the longitudinal ($`T_1`$) and transverse ($`T_2`$) spin relaxation times of conduction band electrons in n-type III-V semiconductors. In particular, we find that the interplay between the Dyakonov-Perel process and an additional spin relaxation channel, which originates from the electron wave vector dependence of the electron $`g`$-factor, yields a maximal $`T_2`$ at a finite magnetic field. We compare our results with existing experimental data on n-type GaAs and make specific additional predictions for the magnetic field dependence of electron spin lifetimes.

The electron spin in a semiconductor is a robust object which can be utilized to add new functionality to existing electronic devices or to even build completely new devices based on this spin degree of freedom. Wolf01  Establishment of successful spintronics devices requires a thorough understanding of the electron spin dynamics in a semiconducting environment. In particular, spin relaxation processes need to be identified and controlled.

Important electron spin relaxation processes in n-type semiconductors include the Elliott-Yafet (EY) process Elliott54 ; Yafet63 , that leads to spin flip scattering and, in semiconductors without inversion symmetry, the Dyakonov-Perel (DP) process Dyakonov72 , in which spin states precess because of spin off-diagonal Hamiltonian matrix elements resulting from a combination of the spin-orbit interaction and inversion asymmetry. Typically, the DP mechanism dominates the spin dynamics in n-type III-V semiconductors. An external magnetic field, in many cases required for control and manipulation in spintronics devices, can significantly influence electron spin dynamics. A magnetic field has two main effects on electron spin relaxation: (i) it quenches the DP process thereby tending to extend the spin lifetimes as a function of the magnetic field Ivchenko73 , and (ii) it opens an additional spin relaxation process which tends to reduce the spin lifetimes in applied magnetic fields. Margulis83  The latter process is due to the wave vector dependence of the conduction band (CB) electron $`g`$-factor. As a result of the variations in the $`g`$-factor, electrons in different quantum states precess about a transverse magnetic field at different rates and thus lose spin coherence. For brevity we will refer to this process as a variable $`g`$-factor (VG) mechanism.

In contrast to previous studies Elliott54 ; Yafet63 ; Dyakonov72 ; Ivchenko73 ; Margulis83 ; Chazalviel75 ; Fishman77 ; Wu00 ; Lau01 ; Song02  of spin relaxation in (bulk) n-type III-V semiconductors, we simultaneously treat the EY, DP, and VG processes on an equal footing and focus on the interplay between the various spin relaxation processes as a function of the magnetic field. Thereby, we are able to study in detail the competition between the quenching of the DP process and the appearance of the VG process.

Specifically, we calculate the longitudinal ($`T_1`$) and transverse ($`T_2`$) spin relaxation times as a function of temperature, electron density and magnetic field. We find that the VG process dominantly influences the transverse ($`T_2`$) spin relaxation time. In particular, as a result of the competition between the quenching of the DP process and the introduction of the VG process, there is a magnetic field for which the transverse ($`T_2`$) spin lifetime is maximal. From the slope of $`T_2`$ at small magnetic fields it is moreover possible to determine whether the DP or the EY process dominates spin relaxation at zero magnetic field. In contrast, the magnetic field dependence of the longitudinal ($`T_1`$) lifetime is essentially unaffected by the VG process and dominated by the quenching of the DP process. Thus, they generally increase with field and saturate at a value given by the EY process.

In an applied magnetic field, the CB electrons in a III-V semiconductor, e.g. GaAs, are described by the Hamiltonian Ogg66 ; Golubev85

$`_{\alpha \beta }(\stackrel{}{K})=ϵ(\stackrel{}{K})\delta _{\alpha \beta }+{\displaystyle \frac{\mathrm{}}{2}}\left[\stackrel{}{\mathrm{\Omega }}_L+\stackrel{}{\mathrm{\Omega }}_{IA}\left(\stackrel{}{K}\right)+\stackrel{}{\mathrm{\Omega }}_g\left(\stackrel{}{K}\right)\right]\stackrel{}{\sigma }_{\alpha \beta },`$ (1)
where $`\stackrel{}{K}=\stackrel{}{k}(e/\mathrm{}c)\stackrel{}{A}(\stackrel{}{r})`$ \[$`\stackrel{}{A}(\stackrel{}{r})`$ is the vector potential\], $`ϵ(\stackrel{}{k})`$ is the Kramers degenerate dispersion of CB electrons, $`\mathrm{}\stackrel{}{\mathrm{\Omega }}_L=\mu _Bg^{}\stackrel{}{B}`$ is the Larmor frequency, $`\mathrm{}\stackrel{}{\mathrm{\Omega }}_{IA}(\stackrel{}{K})=2\delta _0\stackrel{}{\kappa }(\stackrel{}{K})`$ is the splitting of the CB dispersion due to the combination of spin-orbit interaction and inversion asymmetry, and

$`\mathrm{}\stackrel{}{\mathrm{\Omega }}_g(\stackrel{}{K})=2a_4K^2\stackrel{}{B}+2a_5\{\stackrel{}{K},\stackrel{}{B}\stackrel{}{K}\}+2a_6\stackrel{}{\tau }(\stackrel{}{K},\stackrel{}{B})`$ (2)
is a term which gives rise to a wave vector dependence of the CB electron $`g`$-factor. The definitions of the vectors $`\stackrel{}{\kappa }(\stackrel{}{K})`$ and $`\stackrel{}{\tau }(\stackrel{}{K},\stackrel{}{B})`$ and of the parameters $`\delta _0`$, $`a_i,i=4,5,6`$ are given in Refs. Ogg66 ; Golubev85  and $`\{,\}`$ indicates an anticommutator.

Our calculation starts from the full quantum kinetic equations for the contour-ordered Green functions Kinetic , from which we derive, considering a classical homogeneous magnetic field and using the fact that wave vector scattering is essentially instantaneous on the time scale of spin relaxation, a semiclassical kinetic equation for the CB electron density matrix. We then linearize this kinetic equation with respect to the CB electron spin density, assuming, as an initial condition, small spin polarization. Treating scattering processes in the Born approximation and expanding the collision integrals up to second order in the wave vector transfer (diffusion approximation Kinetic ), we finally obtain a generalized Fokker-Planck-Landau equation for the spin density $`\stackrel{}{S}(\stackrel{}{k}t)`$ which, in atomic units (with magnetic field along the z-axis) reads

$`\left[_ti\mathrm{\Omega }_C\widehat{}_z\right]\stackrel{}{S}\left(\stackrel{}{k}t\right)=\left[\stackrel{}{\mathrm{\Omega }}_L+\stackrel{}{\mathrm{\Omega }}_{IA}\left(\stackrel{}{k}\right)+\stackrel{}{\mathrm{\Omega }}_g\left(\stackrel{}{k}\right)\right]\times \stackrel{}{S}\left(\stackrel{}{k}t\right)+\left[\widehat{𝒟}{\displaystyle \frac{1}{2\tau _1\left(k\right)}}\widehat{}^2\right]\stackrel{}{S}\left(\stackrel{}{k}t\right)\mathrm{\Gamma }\left(\stackrel{}{k}\right)\stackrel{}{S}\left(\stackrel{}{k}t\right),`$ (3)
with $`\mathrm{\Omega }_C`$ the cyclotron frequency, $`\widehat{}_z`$ and $`\widehat{}^2`$ the z-component and the squared total angular momentum operator in wave vector space, respectively, $`1/\tau _1(k)`$ the sum of the (on-shell) wave vector relaxation rates for the various scattering processes, and $`\widehat{𝒟}`$ a differential operator in $`k=|\stackrel{}{k}|`$ relevant to inelastic scattering processes.

Equation (3) contains EY, DP, and VG processes and accounts for Larmor precession and orbital motion of the CB electrons in the magnetic field. More specifically, the EY process, due to genuine spin flip scattering events, is given by the tensor $`\mathrm{\Gamma }`$, whereas the DP and VG processes originate from the interplay of spin conserving wave vector scattering events described by the differential operator $`\widehat{𝒟}(1/2\tau _1)\widehat{}^2`$ and the torque forces due to $`\stackrel{}{\mathrm{\Omega }}_{IA}`$ and $`\stackrel{}{\mathrm{\Omega }}_g`$, respectively. The orbital motion encoded in $`i\mathrm{\Omega }_C\widehat{}_z`$ and, to a lesser extend, the torque force due to $`\stackrel{}{\mathrm{\Omega }}_L`$ lead to a quenching of the DP process.

It is possible to derive from Eq. (3) general expressions for the spin relaxation rates without specifying whether the scattering processes are elastic ($`\widehat{𝒟}=0`$) or inelastic ($`\widehat{𝒟}0`$). To that end, we follow Ref. Dyakonov72  and employ a perturbative approach with respect to the torque forces. Our results are therefore valid for $`|\stackrel{}{\mathrm{\Omega }}_{IA}+\stackrel{}{\mathrm{\Omega }}_g|\tau _1<1`$. Expanding $`\stackrel{}{\mathrm{\Omega }}_{IA}(\stackrel{}{k})`$ and $`\stackrel{}{\mathrm{\Omega }}_g(\stackrel{}{k})`$ in terms of spherical harmonics $`Y_{lm}(\mathrm{\Theta },\mathrm{\Phi })`$, we find ($`i=1,2`$Bronold02

$`\left[T_i\right]^1=\left[T_i^{EY}\right]^1+\left[T_i^{DP}\right]^1+\left[T_i^{VG}\right]^1,`$ (4)
with the EY contributions (due to spin-flip scattering)

$`\left[T_1^{EY}\right]^1=2\left[T_2^{EY}\right]^1={\displaystyle \frac{32\pi }{3}}C_{sf}^2{\displaystyle \frac{k^4}{\tau _1(k)}},`$ (5)
the DP contributions (due to inversion asymmetry)

$`\left[T_1^{DP}\right]^1`$ $`=`$ $`4|C_{31}^Y|^2\stackrel{~}{\tau }_{31}^3+4|C_{33}^Y|^2\stackrel{~}{\tau }_{33}^3,`$ (6)
$`\left[T_2^{DP}\right]^1`$ $`=`$ $`2|C_{32}^Z|^2\stackrel{~}{\tau }_{32}^3+2|C_{31}^Y|^2\stackrel{~}{\tau }_{31}^3+2|C_{33}^Y|^2\stackrel{~}{\tau }_{33}^3,`$ (7)
and the VG contributions (due to the wave vector dependence of the CB electron $`g`$-factor)

$`\left[T_1^{VG}\right]^1`$ $`=`$ $`4|D_{21}^Y|^2\stackrel{~}{\tau }_{21}^1,`$ (8)
$`\left[T_2^{VG}\right]^1`$ $`=`$ $`|D_{00}^Z|^2\stackrel{~}{\tau }_{00}^1+|D_{20}^Z|^2\stackrel{~}{\tau }_{20}^1+2|D_{21}^Y|^2\stackrel{~}{\tau }_{21}^3.`$ (9)
Here, $`C_{lm}^i`$ and $`D_{lm}^i`$ ($`i=X,Y,Z`$) are the expansion coefficients of $`\mathrm{\Omega }_{IA}^i(\stackrel{}{k})`$ and $`\mathrm{\Omega }_g^i(\stackrel{}{k})`$, respectively, and ($`\nu =1,2,3`$)

$`\stackrel{~}{\tau }_{lm}^\nu =\text{Re}C_l(k)\tau _{lm}^\nu (k),`$ (10)
where the brackets denote an average over $`k`$ defined as $`(\mathrm{})=_0^{\mathrm{}}𝑑kk^2f(k)\overline{f}(k)(\mathrm{})/4\pi _0^{\mathrm{}}𝑑kk^2f(k)\overline{f}(k)`$, with $`\overline{f}(k)=1f(k)`$ and $`f(k)`$ the equilibrium Fermi distribution function. The generalized wave vector relaxation time $`\tau _{lm}^\nu (k)`$ satisfies a differential equation

$`\left[\widehat{𝒟}+i\left(m\mathrm{\Omega }_C+\mathrm{\Omega }_L^\nu \right){\displaystyle \frac{1}{\tau _l\left(k\right)}}\right]\tau _{lm}^\nu \left(k\right)=C_l\left(k\right),`$ (11)
with $`C_l(k)=C_{IA}k^3\delta _{l3}+C_g(B)k^2[\delta _{l2}+\delta _{l0}]`$, $`\mathrm{\Omega }_L^\nu =\mathrm{\Omega }_L\left[\delta _{\nu 1}\delta _{\nu 2}\right]`$, and $`1/\tau _l(k)=l(l+1)/2\tau _1(k)`$. The constants characterizing the three spin relaxation processes are, respectively, $`C_{sf}=\delta ^2(\mathrm{\Delta }+2ϵ_g)R_0m_0/2\mathrm{\Delta }ϵ_gm^{}`$, $`C_{IA}=2\delta _0/R_0a_0^3`$, and $`C_g(B)=2\mu _BB/R_0a_0^2`$, where $`\delta ^2=2\mathrm{\Delta }^2/(\mathrm{\Delta }+ϵ_g)(2\mathrm{\Delta }+3ϵ_g)`$, $`\mathrm{\Delta }`$ is the spin-orbit splitting, $`ϵ_g`$ is the band gap, $`R_0`$ and $`a_0`$ are the Rydberg energy and the Bohr radius, respectively, $`m^{}`$ and $`m_0`$ are the CB electron mass and the mass of a bare electron, respectively, and $`\mu _B`$ is the Bohr magenton. The detailed form of $`\widehat{𝒟}`$ depends on the scattering processes and does not concern us here. Bronold02

Note, as a consequence of the orthogonality of the angle dependences, the EY, DP, and VG spin relaxation rates are additive. The generalized relaxation rate $`1/\tau _{lm}^\nu (k)`$, on the other hand, is in general not proportional to the sum of the (on-shell) relaxation rates $`1/\tau _l(k)`$ because of inelasticity. A Matthiessen-type rule for $`1/\tau _{lm}^\nu (k)`$ only holds for elastic scattering (see below).

We are interested in the magnetic field dependence of the spin relaxation processes which, at least qualitatively, should not depend on the approximation adopted to decribe the scattering events. In the following, we treat therefore all scattering processes in the elastic approximation and neglect $`\widehat{𝒟}`$ in Eq. (11). Specifically, we take scattering on ionized impurities, acoustic phonons, and longitudinal optical (LO) phonons into account. The elastic approximation restricts our results to low enough temperatures, where electron-impurity scattering dominates, and to high enough temperatures, where electron-phonon scattering becomes essentially elastic.

Within the elastic approximation Eq. (11) reduces to an algebraic equation which is readily solved to yield

$`\stackrel{~}{\tau }_{lm}^\nu ={\displaystyle \frac{\tau _l(k)\left[C_{IA}^2k^6\delta _{l3}+C_g^2(B)k^4(\delta _{l2}+\delta _{l0})\right]}{1+[(m\mathrm{\Omega }_C+\mathrm{\Omega }_L^\nu )\tau _l(k)]^2}}.`$ (12)
The k-average can be obtained either numerically or, at low and high temperatures, with saddle point techniques exploiting the peaked structure of the integrands. Within the elastic approximation it is sufficient to adopt the latter. Details of the calculation will be given elsewhere. Bronold02

In Fig. 1, we show calculated longitudinal ($`T_1`$) and transverse ($`T_2`$) spin relaxation times for GaAs as a function of magnetic field at $`T=0`$ and an electron density of $`n=10^{18}\text{cm}^3`$. We show separately the contributions to the spin relaxation times from the EY, DP and VG processes and the total spin relaxation time including all three spin relaxation processes. In the insets of Fig. 1, we give the total spin relaxation time for various electron densities at $`T=0`$. The parameters needed to specify $`\stackrel{}{\mathrm{\Omega }}_g(\stackrel{}{k})`$ have been previously obtained partly experimentally by measuring combined cyclotron resonances ($`a_4,a_6`$) and partly theoretically within a five-level Kane model ($`a_5`$): $`(a_4,a_5,a_6)=(97,8,49)\times 10^{24}\text{eVcm}^2\text{Oe}^1`$Golubev85  The parameter defining $`\stackrel{}{\mathrm{\Omega }}_{IA}(\stackrel{}{k})`$ is given by $`\delta _0=0.06\mathrm{}^3/\sqrt{(2m^{})^3ϵ_g}`$Aronov83  The remaining parameters, such as the effective CB electron mass or the deformation potential are available from standard data bases. parameter

For the temperature and density conditions in Fig. 1, the electrons are degenerate and electron-ionized impurity scattering dominates. The VG process makes a small contribution to $`T_1`$ which is dominated by the DP process at zero magnetic field. As the magnetic field is increased, the DP process is quenched. Thus, $`T_1`$ increases monotonically with increasing magnetic field saturating at high field at a value determined by the EY process which is not affected by the magnetic field. comment  If the material parameters had been such that the EY process dominated the DP process for $`T_1`$ relaxation at zero magnetic field, $`T_1`$ relaxation would not be significantly affected by the applied field. By contrast the VG process makes a significant contribution to $`T_2`$ relaxation. At small applied magnetic fields the $`T_2`$ lifetime increases with increasing magnetic field, but as the field continues to increase the VG process begins to dominate the relaxation so that $`T_2`$ has a maximum and begins to decrease for larger magnetic fields. If the material parameters had been such that the EY process dominated the DP process for $`T_2`$ relaxation at zero magnetic field, the $`T_2`$ relaxation would monotonically decrease with increasing magnetic field.

The solid squares and triangles in the lower panel of Fig. 1 are measured $`T_2`$ spin lifetimes in GaAs at 5K from Ref. Kikkawa98  at electron densities of $`1\times 10^{18}\text{cm}^3`$ and $`5\times 10^{18}\text{cm}^3`$. (Data for an electron concentration of $`1\times 10^{16}\text{cm}^3`$ at 5K was also presented in Ref. Kikkawa98 , but at this low density the electrons are bound to isolated donors and our theory does not apply.) There is good (order of magnitude) agreement between our calculation and these measured results, although there were no adjustable parameters. Unfortunately, the magnetic fields in Ref. Kikkawa98  are not high enough to capture any effects due to the VG process. In particular, our prediction of the maximum of $`T_2`$ remains to be experimentally verified.

In Fig. 2, we show the various contributions to the $`T_1`$ and $`T_2`$ spin relaxation for GaAs as a function of magnetic field at $`T=100\text{K}`$ and an electron density of $`n=10^{17}\text{cm}^3`$. In the insets of Fig. 2, we show the total spin relaxation time as a function of magnetic field for various temperatures at $`n=10^{17}\text{cm}^3`$. For the temperature and density conditions in Fig. 2, the electrons are non-degenerate and electron-LO-phonon scattering is the dominant scattering process. As for the degenerate electron case, the VG process makes a small contribution to $`T_1`$ which is again dominated by the DP process at zero magnetic field. The DP process is quenched by the field so that $`T_1`$ increases with field at small fields and saturates at a value determined by the EY process at large fields. comment  Similar to the degenerate electron case, the VG process makes a substantial contribution to $`T_2`$ relaxation. At small fields the $`T_2`$ lifetime increases with increasing field and at large fields the VG process begins to dominate the relaxation so that $`T_2`$ has again a maximum at some finite magnetic field. The sign of the slope in $`T_2`$ at small magnetic fields is again a clear signature of whether the EY process ($`T_2`$ decreases with increasing field) or DP process ($`T_2`$ increases with increasing field) dominates $`T_2`$ relaxation at zero magnetic field. Note, the qualitative behavior of the longitudinal and transverse spin relaxation times with increasing magnetic field is similar for degenerate and non-degenerate electrons, but the magnitude of the change is larger for non-degenerate electrons.

In summary, based on a systematic kinetic approach, which treats the EY, DP, and VG processes on an equal footing, we calculated the longitudinal ($`T_1`$) and transverse ($`T_2`$) spin relaxation times of CB electrons in n-type III-V semiconductors as a function of temperature, electron density, and magnetic field. At finite magnetic field, the VG process competes with the DP and EY processes. We find that, as a consequence of the interplay of the DP and the VG processes, $`T_2`$ can have a maximum as a function of magnetic field. In contrast, $`T_1`$ is not affected by the VG process and increases with magnetic field until it saturates at a value determined by the EY process. The sign of the change in $`T_2`$ with increasing magnetic field at small fields indicates, moreover, whether the EY process or the DP process dominates $`T_2`$ relaxation at zero magnetic field. Our calculated results are in good agreement with existing experimental data in n-type GaAs and we make additional specific predictions for the magnetic field dependence of electron spin lifetimes that are subject to experimental check.

This work was supported by the SPINs program of the U.S. Defense Advanced Research Projects Agency.


format: text

text:
# Decoherence-free generation of many-particle entanglement by adiabatic ground-state transitions

## Abstract

We discuss a robust and decoherence insensitive method to create many-particle entanglement in a spin system with controllable collective interactions of the Lipkin-Meshkov-Glick type and propose an implementation in an ion trap. An adiabatic change of parameters allows a transfer from separable to a large variety of entangled eigenstates. We shown that under certain conditions the Hamiltonian possesses a supersymmetry permitting an explicit construction of the ground state at all times of the adiabatic transfer. Of particular interest is a transition in a non-degenerate ground state with a finite energy gap since here the influence of collective as well as individual decoherence mechanisms is substantially reduced. A lower bound for the energy gap is given.

Entanglement is one of the most characteristic features of quantum systems and lies at the heart of quantum information processing and computing entanglement . While in few-particle systems entanglement is by now reasonably well understood entanglement-definition  and practical schemes for its creation and manipulation are developed few-particle-entanglement-experiments ; WinelandScience2000 , many-particle entanglement is still an open field of research with a large unexplored potential for applications. There exist for example proposals to implement quantum computation in an initially entangled many-particle system by performing only measurements and single qubit operations, both of which are relatively easy to implement Briegel . A necessary prerequisite for this is however a specific many-particle entanglement. Since entangled states become increasingly susceptible to environmental interactions if the number of particles increases, an important practical challenge is the design of robust and most importantly decoherence-resistant mechanisms for its generation. We here propose and analyse such a mechanism which is based on adiabatic transitions in a spin system with controllable collective interactions. Robustness against decoherence and short process times are achieved by choosing ground-state transitions with a large energy gap very similar to the ideas used in adiabatic quantum computation Preskill .

Let us consider a collection of $`N`$ interacting spin $`1/2`$ systems described by a generalization of the Lipkin-Meshkov-Glick (LMG) Hamiltonian Lipkin1965

$`H=\xi \left[\lambda \chi _1\chi _2\widehat{J}_z+\chi _1^2\widehat{J}_x^2+\chi _2^2\widehat{J}_y^2+2\mu \chi _2^2\widehat{J}_y\right],`$ (1)
where the $`\widehat{J}_i`$, ($`i\{x,y,z\}`$) are the total-spin operators of the ensemble. We show that this Hamiltonian can be implemented for $`\mu =0`$ with controllable parameters by a generalization of the ion-trap scheme suggested by Sørensen and Mølmer Sorensen1999 , where cold ions interact via a common trap oscillation and are driven by bichromatic lasers fields. In contrast to the effective Hamiltonian of Sorensen1999  the generalized LMG interaction (1) can not be solved exactly. It does provide however the possibility for adiabatic and therefore robust transitions between separable and entangled many-particle states. The most important feature of (1) is that such transitions are possible while staying in a non-degenerate ground state with a finite energy gap and that a large class of entangled target states are accessible by varying $`\mu `$. We show that in the special case of $`\lambda =1`$, the LMG Hamiltonian (1) has a supersymmetry and the ground-state can be explicitly constructed for all times of the transfer process. The gap can be made rather large and hence the process can be fast despite of its adiabatic nature. For the same reason the influence of decoherence is strongly reduced: Collective decoherence due to noise in the external control parameters is eliminated by the adiabatic nature of the process and the influence of independent individual reservoir couplings is suppressed due to the presence of a finite energy gap.

Let us first discuss the case of a negative coupling parameter $`\xi <0`$. If $`\mu =0`$ and $`\lambda N`$, the system described by eq.(1) undergoes a quantum-phase transition Sahdev  when the interaction parameters $`\chi _1`$ and $`\chi _2`$ are changed (case i). This transition can be described analytically in the semiclassical limit with $`J=N/2`$ and $`N1`$. If $`\chi _1=\chi _2`$, $`\widehat{J}_z`$ is a conserved quantity and (1) has a trivial anharmonic spectrum: $`H\xi \chi _2^2\left[\lambda \widehat{J}_z\widehat{J}_z^2+\widehat{\text{J}}^2\right].`$ The ground state of the system has maximum total angular momentum $`J=N/2`$ and a $`z`$ projection $`m_z=\pm N/2`$ depending on the sign of $`\lambda `$. Both of these states denoted by $`|\mathrm{}`$ and $`|\mathrm{}`$ are separable many-particle states. On the other hand in the limit $`\chi _1=0`$ the terms in (1) containing $`\widehat{J}_x`$ and $`\widehat{J}_z`$ vanish and the Hamiltonian approaches $`H\xi \chi _2^2\widehat{J}_y^2.`$ It has again a trivial spectrum, however with two degenerate ground states $`|m_y=\pm N/2`$. Both of these, taken individually, are separable. Symmetric or anti-symmetric superpositions of them form however the $`N`$-spin analogues of the Greenberger-Horne-Zeilinger states (GHZ states). Using the Schwinger-representation of angular momenta Schwinger  one finds that the state $`|m_z=N/2`$ is adiabatically connected only to one particular superposition due to the symmetry of the Lipkin interaction. Thus

$$|\mathrm{}\stackrel{(i)}{}\frac{1}{\sqrt{2}}[|++\mathrm{}+\mathrm{e}^{i\pi J}|\mathrm{}],$$
(2)
which corresponds to the generation of the $`N`$-particle analog of the GHZ state in the $`\sigma _y`$ basis. $`|\pm `$ denote single-particle eigenstates of $`\sigma _y`$ with $`m_y=\pm 1/2`$ respectively. Due to the phase factor $`e^{i\pi J}`$ the entangled state depends sensitively on the total number of atoms even in the limit $`N\mathrm{}`$ Pokrovsky .

For a positive coupling parameter, $`\xi >0`$, three cases need to be distinguished: $`\mu =0`$ and the total number of spins is odd (case ii); $`\mu =0`$ and the total number of spins is even (case iii); and $`\mu 0`$ (case iv). In all cases the ratio $`\chi _1/\chi _2`$ is again rotated from unity to zero.

Since $`\xi >0`$ the initial ground state in cases (ii) and (iii) is $`|m_z=\pm N/2`$ depending on the sign of $`\lambda 0`$ and is separable. The final Hamiltonian is again $`H=\xi \chi _2^2\widehat{J}_y^2`$, whose ground state is now however the eigenstate of $`\widehat{J}_y`$ with smallest value of $`|m_y|`$. Thus for an odd number of particles (case ii) there are two degenerate ground states with $`m_y=\pm 1/2`$, both of them being maximally entangled. Making use of the symmetry of the interaction and using the Schwinger-representation we find that the adiabatic transition leads to the mapping

$$|\mathrm{}\stackrel{(ii)}{}\frac{1}{\sqrt{2}}[|(+)^n()^{n1}_\mathrm{s}+i|(+)^{n1}()^n_\mathrm{s}],$$
(3)
where $`n=N/2`$ and the subscript “s” denotes symmetrization.

In case (iii) the final ground state ($`\chi _10`$) is non-degenerate, has spin projection $`m_y=0`$ and is maximally entangled.

$$|\mathrm{}\stackrel{(iii)}{}|m_y=0=|(+)^{N/2}()^{N/2}_\mathrm{s}.$$
(4)
In this case there is no merging of eigenstates and consequently no phase transition. The absence of degeneracy during the entire adiabatic transfer makes case (iii) particularly interesting because here decoherence is strongly suppressed by the presence of a finite energy gap. The target state is however fixed to the special case $`|m_y=0`$.

The variety of accessible target states in a non-degenerate adiabatic ground-state transition can be substantially increased by adding a linear interaction proportional to $`\widehat{J}_y`$ to the LMG Hamiltonian, i.e. by allowing for a non-vanishing value of $`\mu `$ in eq.(1), (case iv): If we assume for simplicity that $`\mu =m`$, with $`m`$ being an integer or half-integer with $`m\{J,J+1,\mathrm{},J\}`$, the final ($`\chi _10`$) Hamiltonian approaches $`H\xi \chi _2^2\left[\widehat{J}_y^2+2m\widehat{J}_y\right]`$. Its ground state is $`|m_y=m`$ and thus can be adjusted to have any eigenvalue of $`\widehat{J}_y`$, which is maximally entangled unless $`m=\pm J`$.

$$|\mathrm{}\stackrel{(iv)}{}|m_y=m.$$
(5)
In all four cases considered an adiabatic change of $`\chi _1/\chi _2`$ from unity to zero leads to the generation of a well defined entangled state. In cases (iii) and (iv) this is moreover possible with a finite energy gap. The adiabatic process is neither of the Landau-Zener nor the STIRAP type VitanovAdvAMO2001 . We have illustrated the different scenarios for even $`N`$, $`\mu =0`$, and a positive linear term in (1) in Fig. 1.

To find the ground-state of (1) for arbitrary values of $`\chi _1/\chi _2`$ as e.g. during the transfer process is a difficult task. We will show now that (1) possesses a super-symmetry (SUSY) for $`\lambda =1`$, allowing for an explicit construction of the ground state for cases (iii) and (iv). (The existence of an extra symmetry has been suspected GargPRB2001  but has not been understood so far.) For $`\lambda =1`$ the Hamiltonian can be factorized as

$`{\displaystyle \frac{H}{\xi }}`$ $`=\left(\chi _1\widehat{J}_x+i\chi _2\widehat{J}_yi\chi _2\mu \right)\left(\chi _1\widehat{J}_xi\chi _2\widehat{J}_y+i\chi _2\mu \right)`$ (6)
$`\chi _2^2\mu ^2.`$
Due to the SUSY the spectrum consists of twofold degenerate and a nondegenerate state WittenNuclPhys1981 . For $`\mu =m\{J,\mathrm{},J\}`$ the ground state is nondegenerate and obeys $`(\chi _1\widehat{J}_xi\chi _2\widehat{J}_y+i\chi _2m)|\psi _0=0`$. One finds

$`|\psi _0=𝒩\mathrm{exp}\left(\gamma \widehat{J}_z\right)|m_y=m`$ (7)
where $`𝒩`$ is a normalization constant and $`\mathrm{tanh}(\gamma )=\chi _1/\chi _2`$ for $`\chi _1\chi _2`$. $`|\psi _0`$ is an entangled state of the $`N`$ spins for any value of $`\gamma 0`$. I.e. by varying $`\gamma `$ and choosing $`m`$, we have access to a rich variety of entangled many-particle states while staying in the lowest energy state.

We proceed by discussing the conditions for adiabaticity and the sensitivity of the process to decoherence. In this context particular attention has to be given to the phase transitions in cases (i) and (ii) since they are associated with a merging of pairs of energies. Due to the symmetry of (1) only one superposition of the associated states is however coupled to the nondegenerate initial ground state. Consequently for the transition to be adiabatic it is sufficient that the characteristic time of the transfer $`T`$ is much larger than the typical inverse frequency difference $`\mathrm{}/\mathrm{\Delta }E`$ to the next excited states. For the same reason collective decoherence processes caused e.g. by fluctuations in the external parameters $`\xi ,\lambda ,\chi _1`$, and $`\chi _2`$ are suppressed in the present system by an exponential Boltzmann factor $`\mathrm{exp}(\beta \mathrm{\Delta }E)`$, where $`\beta ^1=k_BT`$ is the thermal energy of the heat bath. Decoherence processes caused by independent heat bath couplings of the individual spins do not have the symmetry of (1) however and hence do couple the degenerate states in cases (i) and (ii). Only in cases (iii) and (iv), i.e. without a quantum phase transition, there is always a finite energy gap to all other states with the same total angular momentum $`J`$. Thus the entanglement generation in cases (iii) and (iv) will be robust against collective and individual decoherence processes, provided the energy gap is sufficiently large.

It is not possible to give an analytic expression for the energy gap $`\mathrm{\Delta }E`$ in the general case. Numerical investigations for up to 50 particles indicate that for a transfer efficiency close to 100% it is sufficient that

$$\chi _1^{\mathrm{max}}\sqrt{T}\chi _2^{\mathrm{max}}\sqrt{T}1$$
(8)
with $`\lambda N`$ for case ($`i`$) and $`\lambda 1`$ for cases ($`iiiv`$). $`T`$ is the characteristic transfer time. An estimate for $`\mathrm{\Delta }E`$ in cases (ii) and (iii), i.e. for $`\mu =0`$ can be obtained as follows: Using a variational method with a trial functions $`|\mathrm{\Psi }^{(N)}`$ one finds an energy estimate $`H_N`$ for the ground state with $`J=N/2`$, $`E_0^{(N)}H_N`$. Secondly one can apply Temple’s formula Temple52  to obtain a lower bound $`E_0^{(N)}H_N\frac{\mathrm{\Delta }H^2_N}{E_1^{(N)}H_N}`$, where $`\mathrm{\Delta }H^2_N=\left(HH\right)^2_N`$ is the energy fluctuation in the trial state, $`E_1^{(N)}`$ is the energy of the first excited state. One can show that Temple’s formula gives the best lower bound when using only $`\mathrm{\Delta }H^2_N`$, $`H_N`$ and $`E_1^{(N)}`$ as parameters. Furthermore we make use of the inequality $`E_1^{(N)}E_0^{(N2)}`$ between the energy of the first excited state for $`N`$ particles and the ground state for $`N2`$ particles, in both cases with maximum $`J`$, which can easily be proven. Applying this inequality iteratively leads to

$$\frac{\mathrm{\Delta }E^{(N)}}{H_{N2}H_N}\frac{1}{2}\left(1+\sqrt{14A}\right)$$
(9)
where

$$A>\frac{\mathrm{\Delta }H^2_N}{\left(H_{N2}H_N\right)\left(H_{N4}H_{N2}\right)}.$$
(10)
By choosing the trial function of the ground state close to the exact one it is possible to achieve $`4A1`$, i.e. $`\mathrm{\Delta }E^{(N)}H_{N2}H_N`$. For the simple trial function $`|\mathrm{\Psi }^{(N)}=\alpha _1|m_y=0+\alpha _2|m_z=N/2`$ one finds $`H_{N2}H_N=\beta \lambda \left|\xi \right|\chi _1\chi _2|_{\mathrm{max}}`$, where $`\beta `$ is a numerical factor of order unity, which varies only very slowly with $`N`$. Thus a reasonable estimate for the energy gap in case ($`iii`$) is given by $`\beta \lambda \left|\xi \right|\chi _2^2`$. This is also confirmed by our numerical calculations for particle numbers up to $`N=50`$. If $`\beta \lambda \left|\xi \right|\chi _2^2`$ is sufficiently larger than the thermal energy of the environment, the probability of decoherence processes is strongly suppressed.

Let us now discuss a possible implementation of the Lipkin Hamiltonian (1) with $`\mu =0`$ in an ion-trap system. Consider a linear trap with a string of ions with two relevant internal levels $`|g`$ and $`|e`$. The ions are assumed to be cooled such that only the in-phase collective oscillation of all ions is excited. The corresponding oscillation frequency is denoted by $`\nu `$. The two internal levels of the ions are coupled by two laser fields with frequencies $`\omega _1`$ and $`\omega _2`$, and slowly-varying Rabi-frequencies $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$. Assuming that both fields couple all ions in the same way, we can describe the system by the Hamiltonian $`H=H_0+H_{\mathrm{int}}`$, where $`H_0=\mathrm{}\nu \widehat{c}^{}\widehat{c}+\mathrm{}\omega _{eg}\widehat{J}_z`$, $`\widehat{c}`$ and $`\widehat{c}^{}`$ being the annihilation and creation operators of the trap oscillation. $`\mathrm{}\omega _{eg}`$ is the energy separation between the two internal states and $`\widehat{J}_z=\frac{1}{2}_{i=1}^N(\sigma _{ee}^i\sigma _{gg}^i)`$, with $`\sigma _{\mu \mu }^i=|\mu _{ii}\mu |`$ being the projector to the internal state $`|\mu `$ of the $`i`$th ion. The interaction Hamiltonian $`H_{\mathrm{int}}`$ is given in rotating wave approximation by

$$H_{\mathrm{int}}=\widehat{J}_+\mathrm{e}^{i\eta (\widehat{c}+\widehat{c}^{})}\left[\mathrm{\Omega }_1\mathrm{e}^{i\omega _1t}+\mathrm{\Omega }_2\mathrm{e}^{i\omega _2t}\right]+h.c.$$
(11)
where $`\eta `$ is the Lamb-Dicke parameter. We assume that the laser frequencies have equal and opposite detuning $`\delta `$ from resonance $`\omega _{1,2}=\omega _{eg}\pm \delta `$. $`\delta `$ is assumed to be large compared to the linewidth of the resonance but sufficiently different from the frequency of the trap oscillation, i.e. $`|\delta |,|\nu \pm \delta |\gamma `$. As a consequence the dominant processes are two-photon transitions leading to a simultaneous excitation of pairs of ions as indicated in Figure 2. For $`\mathrm{\Omega }_1=\mathrm{\Omega }_2`$ this scheme has first been considered by Sørensen and Mølmer Sorensen1999 ; Molmer1999  in the context of quantum computation and dynamical entanglement generation.

We now assume that the ion trap is in the Lamb-Dicke limit, i.e. that the ions are cooled sufficiently enough, such that for all relevant excitation numbers $`n`$ of the trap oscillation $`(n+1)\eta ^21`$ holds. In this limit one can expand the exponent in (11) to first order in $`\eta `$. For large values of $`|\delta |`$ it is convenient to consider this interaction in terms of a coarse-grained Hamiltonian which neglects the effects of rapidly oscillating terms. Using the time-averaging method of ref. James  one arrives at an effective Lipkin Hamiltonian (1) with the identifications $`\xi =(2\nu \eta ^2)/(\delta ^2\nu ^2)`$, $`\lambda =2/\xi \delta `$, and $`\chi _{1,2}=\mathrm{\Omega }_1\mathrm{\Omega }_2`$. $`\mathrm{\Omega }_1=\mathrm{\Omega }_2`$ corresponds to the case $`\chi _1=\mu =0`$ in (1) which is exactly solvable and has been discussed in Refs.Sorensen1999 ; Molmer1999 . It has also been shown in Sorensen1999 ; Molmer1999  that the effective Hamiltonian describes correctly the dynamics of the ions in a coarse-grained time picture. Furthermore the coupling scheme has successfully been implemented in experiments to generate entanglement between 4 ions WinelandScience2000 .

In Fig.3 we have shown as an example the effective dynamics of a system of 4 ions driven by fields $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ for the cases (i) and (iii) following eq.(1) with $`\mu =0`$. One recognizes that a nearly perfect transfer is possible.

In summary we have shown that it is possible to generate specific entangled many-particle states in an ensemble of spins interacting through a collective coupling of the Lipkin-Meshkov-Glick type by adiabatic ground-state transitions. Two scenarios (i,ii) involve a two-fold degenerate ground state at some stages while two others (iii) and (iv) always have a nondegenerate ground state. In all cases there is a finite energy gap to other excited states. This gap can be rather large and thus fast processes are possible despite the required adiabaticity. In the asymptotic limits of the adiabatic transfer, the spectrum and eigenstates of the LMG-Hamiltonian can be exactly calculated. Furthermore for $`\lambda =1`$ there is a supersymmetry allowing for an explicit construction of the ground state for all times. Due to adiabaticity, all transitions are robust against parameter variations. Furthermore due to the symmetry of the coupling and the finite energy gap from the (degenerate or nondegenerate) ground state to excited states, collective decoherence processes are suppressed. In addition in the case of a nondegenerate ground state (case iii and iv) also individual decoherence processes are suppressed. The collective spin Hamiltonian can be implemented in a cold ensemble of ions in a linear trap driven by a nearly-resonant bichromatic fields.

The authors would like to thank C. Mewes and L. Plimak for their help with numerics. This work has been supported by the Deutsche Forschungsgemeinschaft under contract no. Fl-210/10. RU also acknowledges financial support by the Alexander-von-Humboldt foundation.


text:
# Enhanced gauge symmetry from ‘toric’ 𝐺₂ cones

## 1 Introduction

M-theory suggests phenomenologically interesting constructions based on compactification from eleven dimensions. Among these, compactification on 7-dimensional spaces of $`G_2`$ holonomy deserves special attention, since it leads naturally to $`N=1`$ supersymmetric field theories in four dimensions.

It is well-known, however, that compactification on smooth $`G_2`$ spaces has major phenomenological drawbacks, as it fails to produce nonabelian gauge groups and chiral matter . As pointed out in , this can be overcome by considering singular $`G_2`$ backgrounds. In such cases, nonabelian gauge symmetry arises through $`M2`$-branes wrapping some vanishing two-cycles , while net chirality can occur for purely topological reasons . This mechanism was illustrated in  by considering a specific class of conical $`G_2`$ spaces, some of which admit a IIA description in terms of intersecting D6-branes. In particular models, the appearance of chiral matter can also be inferred directly from IIA constructions . Such type II systems consist of D6-branes/orientifold 6-planes, thus admitting a lift to an $`M`$-theory background of $`G_2`$ holonomy.

The class of spaces studied in  is obtained as follows . Starting with an Einstein self-dual space $`M`$ (of positive scalar curvature), one considers its six-dimensional twistor space $`Y`$, which can be written as the sphere bundle associated with the vector bundle $`\mathrm{\Lambda }^{2,}T^{}M`$ of antiselfdual two forms. The twistor space carries a Kahler-Einstein metric:

$$d\rho ^2=d\sigma ^2+|d_A\stackrel{}{u}|^2,$$
(1)
where $`d\sigma ^2`$ is the metric on $`M`$, $`A`$ is the connection induced on $`\mathrm{\Lambda }^{2,}T^{}M`$ by the Levi-Civita connection of $`M`$, and $`\stackrel{}{u}=(u^1,u^2,u^3)`$ are coordinates on the $`S^2`$ fiber, subject to the constraint $`|\stackrel{}{u}|^2=1`$. To obtain a $`G_2`$ space, one takes the metric cone over $`Y`$, where the latter is endowed with the modified metric $`d\rho ^2=\frac{1}{2}(d\sigma ^2+\frac{1}{2}|d_A\stackrel{}{u}|^2)`$. In fact, this cone admits a one-parameter family of $`G_2`$ deformations given by :

$$ds_{G_2}^2=\frac{1}{1(r_0/r)^4}dr^2+\frac{r^2}{2}(d\sigma ^2+\frac{1}{2}(1(r_0/r)^4)|d_A\stackrel{}{u}|^2),r_00.$$
(2)
The conical limit is reached for $`r_0=0`$. It is clear from this expression that the singularities of the conical $`G_2`$ metric are determined by the singularities of $`Y`$: except for the apex, the former are obtained by taking the cone over the latter (and their singularity type coincides with that in $`Y`$). The importance of studying such singularities is clear since they determine the essential physics of M-theory on our backgrounds .

To specify the $`G_2`$ space completely, one must know the (positive curvature) Einstein-selfdual metric on $`M`$. If $`M`$ is smooth, then there are only two choices, namely $`S^4`$ and $`\mathrm{𝐂𝐏}^2`$. For these cases, M-theory physics on the associated $`G_2`$ spaces was studied in . Allowing for singularities in $`M`$ leads to many other examples. Among these are the models considered in , which correspond to $`M=W\mathrm{𝐂𝐏}_{p,q,r}^2`$, endowed with the Einstein-selfdual metric obtained implicitly in . In these examples, the singularities of the twistor space can be determined by elementary methods . However, this becomes a rather difficult task if one wishes to consider more general situations.

A characteristic feature of the examples of  is the presence of a two-torus of isometries. It is by now well-known  that there exist infinitely many inequivalent compact ESD orbifolds111By this we mean a local orbifold (V-manifold). The spaces under consideration are generally not global quotients of some manifold by a finite group. admitting two commuting isometries. Such spaces provide a vast generalization of the models discussed in . In , we developed a method for analyzing the singularities of the twistor space (and thus of the associated $`G_2`$ cone and of its deformation (2)) for this much larger class of ‘toric’ ESD orbifolds.

The approach of  relies on a correspondence between ESD orbifolds and certain hyperkahler cones (a hyperkahler cone is a hyperkahler space which is the metric cone over a compact Riemannian space). Recall that the ESD property of the metric amounts to the quaternion-Kahler condition in four dimensions. Now, it is well-known that every quaternion-Kahler space $`M`$ has an associated hyperkahler cone222The correspondence of  holds irrespective of the sign of the scalar curvature of $`M`$. In the case of interest for us (namely when $`M`$ has positive scalar curvature) the hyperkahler metric on $`X`$ will be Riemannian (positive-definite). Negative curvature quaternion-Kahler spaces have pseudo-Riemannian hyperkahler cones, a situation which is of interest in supergravity. $`X`$ , and this correspondence can be ‘inverted’ i.e. $`M`$ can be re-constructed given $`X`$. This allows one to translate problems in quaternion-Kahler geometry into the language of hyperkahler spaces, which tends to be more amenable to a solution. In our case, the hyperkahler cone is real eight-dimensional and will admit a $`T^2`$’s worth of isometries, induced from the isometries of $`M`$. $`4n`$-dimensional hyperkahler spaces with an $`n`$-torus of isometries have been considered in the mathematics literature  and are called toric hyperkahler. Our particular situation fits into that theory for $`n=2`$. Such spaces can be described as hyperkahler quotients by certain torus actions, much in the spirit of usual toric geometry (which is concerned with Kahler toral quotients). This allows one to reduce geometric questions for $`X`$ to problems in integral linear algebra and combinatorial convex geometry.

The twistor space of $`M`$ can be obtained from $`X`$ by performing a certain $`U(1)`$ Kahler quotient at a positive moment map level. This description allows one to extract the singularities of $`Y`$ by first determining the singularities of $`X`$ and then studying the effect of this Kahler reduction . The outcome of this analysis is a simple algorithm for identifying the singularities of the twistor space, which we summarize in the next section. In  we also used geometric arguments to describe the type IIA reduction along one of the $`U(1)`$ isometries of the associated $`G_2`$ space. For a generic choice of the reduction isometry, the resulting IIA background is strongly coupled, and does not seem to admit a simple perturbative description. However, certain models admit a ‘good isometry’, which leads to a system of $`D6`$-branes upon reduction. Such an isometry (when it exists) can be found by a simple criterion discussed in . Since the IIA metric inherits a $`U(1)`$ isometry induced from the M-theory solution, one can T-dualize along its direction to obtain a IIB background. This will correspond to a system of delocalized 5-branes if a ‘good isometry’ is used for IIA reduction. While models with a good isometry admit a simple IIA/IIB description, we stress that there is no obvious reason to discard other models. From the M-theory perspective, models without a good isometry are as good as any other.

The abstract arguments of  were confirmed in  by direct calculation. Using the result of , we wrote down the $`G_2`$ metric obtained from a general ESD space of positive scalar curvature and admitting a $`T^2`$ of isometries. By performing the reduction, we found the corresponding T-dual IIA and IIB backgrounds. Guided by the analysis of  and the topological arguments of , we computed the explicit asymptotics of relevant fields (metric, coupling constant, and R-R/NS-NS forms) in the vicinity of the branes and determined the relevant RR fluxes. The main results of  and  are briefly reviewed below.

## 2 Singularities of the twistor space

The toric hyperkahler cones associated with our ESD spaces are obtained as hyperkahler quotients of some affine quaternion space $`𝐇^n`$ by an appropriate torus action. Namely, we have $`X=𝐇^n///_0U(1)^{n2}`$, where the subscript $`0`$ indicates that the hyperkahler quotient is taken at zero moment map levels. The action of $`U(1)^{n2}`$ on the quaternion coordinates $`u_1\mathrm{}u_n𝐇`$ is given by:

$$u_k\underset{\alpha =1}{\overset{n2}{}}\lambda _\alpha ^{q_k^{(\alpha )}}u_k,$$
(3)
where $`\lambda _\alpha `$ are complex numbers of unit modulus. The $`(n2)\times n`$ matrix333One has to impose some mild conditions on $`Q`$ to guarantee that the $`U(1)^{n2}`$ action on $`𝐇^n`$ is effective. See  and footnote 6 of this letter. $`Q_{\alpha k}=q_k^{(\alpha )}`$ determines the cone $`X`$. Following toric geometry procedure, we introduce a $`2\times n`$ matrix $`G`$ whose rows form an integral basis for the kernel of $`Q`$. The columns of $`G`$ are two-dimensional integral vectors $`\nu _1\mathrm{}\nu _n`$, which we call toric hyperkahler generators. Splitting every quaternion coordinate into its complex components:

$$u_k=w_k^{(+)}+𝐣w_k^{()},$$
(4)
(where $`𝐢,𝐣,𝐤`$ are the imaginary quaternion units) we obtain the following re-writing of the torus action:

$$w_k^{(+)}\underset{\alpha =1}{\overset{n2}{}}\lambda _\alpha ^{q_k^{(\alpha )}}w_k^{(+)},w_k^{()}\underset{\alpha =1}{\overset{n2}{}}\lambda _\alpha ^{q_k^{(\alpha )}}w_k^{()}.$$
(5)
The twistor space is obtained as a $`U(1)`$ Kahler quotient $`Y=X//_\zeta U(1)`$ of $`X`$ at a positive level $`\zeta `$ (the precise choice of $`\zeta `$ fixes the overall scale of $`Y`$ and $`M`$). The $`U(1)`$ action is given by:

$$w_k^{(\pm )}\lambda w_k^{(\pm )}$$
(6)
and will be called ‘the projectivizing action’. For simplicity, we shall assume 444If this assumption does not hold, then the projectivizing $`U(1)`$ has a trivially acting $`𝐙_2`$ subgroup. In this case, there are two toric hyperkahler cones associated with $`Y`$ and the results described below must be slightly modified. The required modifications are described in . that this action is effective on $`X`$.

In , we show that the singular locus of $`Y`$ is a subset of the so-called distinguished locus $`Y_D`$. The latter consists of two kinds of holomorphically embedded two-spheres, which are distinguished by their position with respect to the $`S^2`$ fibration $`YM`$:

1) $`Y_D`$ contains $`n`$ vertical spheres $`Y_j`$ ($`j=1\mathrm{}n`$), which are fibers of $`Y`$ over $`M`$.

2) it also contains $`2n`$ horizontal spheres $`Y_e`$, which are lifts of spheres lying in $`M`$.

The meaning of the index $`e`$ will become apparent in a moment (it corresponds to the edges of a certain polygon). The union $`Y_V`$ of all vertical spheres and the union $`Y_H`$ of horizontal spheres will be called the vertical and horizontal loci.

A basic observation (which goes back to , see also ) is that the hyperkahler moment map constraints are solved by points $`u𝐇^n`$ whose complex coordinates are solutions of the system:

$$\frac{1}{2}(|w_k^{(+)}|^2|w_k^{()}|^2)=\nu _ka,w_k^{(+)}w_k^{()}=\nu _kb,\mathrm{for}\mathrm{all}k=1\mathrm{}n,$$
(7)
for some parameters $`a=(x_1,y_1)𝐑^2`$ and $`b=(x_3+ix_2,y_3+iy_2)𝐂^2`$. Using this, one shows  that the distinguished locus can be described in terms of the characteristic polygon:

$$\mathrm{\Delta }=\{a𝐑^2|\underset{k=1}{\overset{n}{}}|\nu _ka|=\zeta \}.$$
(8)
This is a planar convex polygon with $`2n`$ vertices, which is symmetric under the point reflection $`aa`$. In particular, $`\mathrm{\Delta }`$ has $`n`$ principal diagonals $`D_j=\{a𝐑^2|a\nu _j=0\}`$. These are the diagonals which pass through the origin and thus connect opposite vertices.

Each vertical sphere $`Y_j`$ is defined by the condition $`u_j=0`$ (i.e. $`w_j^{(+)}=w_j^{()}=0`$) with $`j=1\mathrm{}n`$. The horizontal spheres $`Y_e`$ are indexed by the edges $`e`$ of $`\mathrm{\Delta }`$ and are determined as follows. Given such an edge, define the signs $`ϵ_j(e)=sign(\nu _jp_e)`$, where $`p_e`$ is the two-vector from the origin of the plane to the middle point of $`e`$. Then $`Y_e`$ is the locus in $`Y`$ given by the equations $`w_1^{(ϵ_1)}=w_2^{(ϵ_2)}=\mathrm{}=w_n^{(ϵ_n)}=0`$. As discussed in , $`Y_j`$ is a circle fibration over the principal diagonal $`D_j`$ (the circle fibers collapse to points above the two vertices connected by the diagonal). Each horizontal sphere $`Y_e`$ is a circle fibration over the edge $`e`$, whose circle fibers collapse to points above the vertices connected by this edge. The twistor space admits an antiholomorphic involution acting along its $`S^2`$ fibers. Its restriction to $`Y_H`$ intertwines the horizontal spheres $`Y_e`$ and $`Y_e`$, while covering the point reflection $`\iota :aa`$ of $`\mathrm{\Delta }`$. Finally, the two horizontal spheres associated with adjacent edges of $`\mathrm{\Delta }`$ and the vertical sphere associated with the principal diagonal passing through their common vertex have precisely one point in common. Such special points of the distinguished locus are in one-one correspondence with the vertices of $`\mathrm{\Delta }`$. This situation is summarized in figure 1.

As mentioned above, the singular locus of the twistor space is a subset of its distinguished locus. Namely, some of the vertical and horizontal spheres will consist of cyclic singularities, while some of the special points will carry singularities of their own. The relevant singularity types are described by the following rules:

(a) Given a horizontal sphere $`Y_e`$, consider the integral vector 555With our assumptions on $`Q`$, this vector can never vanish.:

$$\nu _e=\underset{k=1}{\overset{n}{}}ϵ_k(e)\nu _k,$$
(9)
where the signs $`\{ϵ_k(e)\}`$ are defined as explained above. Then $`Y`$ has a $`\mathrm{\Gamma }_e=𝐙_{m_e}`$ quotient singularity along $`Y_e`$, where $`m_e`$ is the greatest common divisor of the two entries of the vector $`\nu _e`$. The orbifold action on the coordinates transverse to $`Y_e`$ can be determined as explained in .

(b) Given a vertical sphere $`Y_j`$, consider the matrix $`\stackrel{~}{Q}_j`$ obtained by deleting the $`j^{th}`$ and $`(j+n)^{th}`$ columns of the $`(n1)\times (2n)`$ matrix $`\stackrel{~}{Q}=\left[\begin{array}{ccc}Q& & Q\\ 1& \mathrm{}& 1\end{array}\right]`$. Then the singularity group $`\mathrm{\Gamma }_j`$ of $`Y`$ along $`Y_j`$ coincides with $`𝐙_{m_j}`$ or $`𝐙_{2m_j}`$, where $`m_j`$ is the greatest common divisor of the components of $`\nu _j`$. To find which of these cases occurs, one computes the integral Smith form666Given an integral $`r\times n`$ matrix $`F`$ with $`rn`$, one can find matrices $`UGL(r,𝐙)`$ and $`VGL(n,𝐙)`$ such that the matrix $`F^{ismith}=U^1FV`$ (the integral Smith form of $`F`$) has the form $`[D,0]`$. Here $`D=diag(t_1\mathrm{}t_r)`$ with $`t_1\mathrm{}t_r`$ some non-negative integers satisfying the division relations $`t_1|t_2|\mathrm{}|t_r`$. These integers are called the invariant factors of $`F`$. The mild condition mentioned in footnote 3 is that all invariant factors of $`Q`$ be equal to $`1`$. of the matrix $`\stackrel{~}{Q}_j`$:

$$\stackrel{~}{Q}_j^{ismith}=[diag(1\mathrm{}1,t_j),0],$$
(10)
where $`t_j=m_j`$ or $`t_j=2m_j`$. The singularity group $`\mathrm{\Gamma }_j`$ coincides with $`𝐙_{t_j}`$. The orbifold action on the transverse coordinates can be determined as explained in .

Note that some of the vertical and horizontal spheres will be smooth (this happens respectively when $`m_e=1`$ and $`t_j=1`$). The singularity type at the special points (which correspond to the vertices of $`\mathrm{\Delta }`$) can also be determined by simple criteria, which the interested reader can find in .

The criteria listed above determine the gauge groups produced at the singular loci of the $`G_2`$ cone: a $`𝐙_m`$ singularity corresponds to an $`SU(m)`$ gauge factor. On the other hand, the second Betti number of our $`G_2`$ cone equals $`n1`$. As in , M-theory on such cones produces an Abelian $`U(1)^{n1}`$ factor, and chiral fermions charged under this group and localized at the cones’ apex. The second Betti number reduces to $`n2`$ (and the Abelian gauge factor is Higgsed down to $`U(1)^{n2}`$) when deforming the cone to the solution with non-vanishing $`r_0`$. All of this is qualitatively identical with the behavior of the models discussed in .

Dimensionally reducing M-theory through one of the $`T^2`$ isometries, one obtains a type IIA background with a non-vanishing RR one-form. Such an isometry is called good if its fixed point set coincides with the singular locus of the $`G_2`$ cone. If a good isometry exists, then the associated type IIA reduction gives a system of D6-branes. It is shown in  that a good isometry exists if all two-vectors $`\nu _e,\nu _j`$ are collinear; in this case, the good isometry corresponds to their common direction, viewed as a direction in the Lie algebra of the two-torus of isometries of the $`G_2`$ space.

Since the resulting IIA background inherits a $`U(1)`$ isometry, one can T-dualize to a IIB solution. It turns out that the orbits of the T-dualizing isometry lie along the $`S^1`$ fibers of the loci $`Y_jD_j`$ and $`Y_ee`$, both of which descend to the IIA solution. In the good isometry case, it follows that T-duality acts along the worldvolume directions, thereby producing delocalized 5-branes in type IIB. This abstract argument can be substantiated by a direct analysis of the relevant field configurations, which was carried out in  and is briefly reviewed below.

## 3 $`G_2`$ metrics and reduction to type II string theory

As discussed above, the ESD space $`M`$ and its associated $`G_2`$ cone, as well as the hyperkahler cone $`X`$ admit a two-torus of isometries for our class of models. The most general ESD metric admitting two commuting isometries was given explicitly in recent work of Calderbank and Pedersen . To describe the solution, one introduces coordinates $`\varphi `$, $`\psi `$ along the $`T^2`$ fibers of $`M`$ and $`\rho ,\eta `$ along the base of its $`T^2`$ fibration. There exists an $`SL(2,𝐙)`$ freedom in the choice of $`\varphi `$ and $`\psi `$, which is related to modular transformations of the lattice $`𝐙^2`$ containing the toric hyperkahler generators. As explained in , the base coordinates can be chosen such that $`\rho 0`$ and $`\eta 𝐑`$, and parameterize the upper half plane (the hyperbolic plane $`^2`$). The metric of  takes the form:

$`d\sigma ^2={\displaystyle \frac{F^24\rho ^2\left(F_\rho ^2+F_\eta ^2\right)}{4F^2}}{\displaystyle \frac{d\rho ^2+d\eta ^2}{\rho ^2}}+{\displaystyle \frac{\left[\left(F2\rho F_\rho \right)\alpha 2\rho F_\eta \beta \right]^2+\left[2\rho F_\eta \alpha +\left(F+2\rho F_\rho \right)\beta \right]^2}{F^2\left[F^24\rho ^2\left(F_\rho ^2+F_\eta ^2\right)\right]}},`$ (11)
where

$$F=\underset{k=1}{\overset{n}{}}\frac{\sqrt{(\nu _k^2)^2\rho ^2+(\nu _k^2\eta +\nu _k^1)^2}}{\sqrt{\rho }}$$
(12)
and $`\alpha =\sqrt{\rho }d\varphi `$, $`\beta =(d\psi +\eta d\varphi )/\sqrt{\rho }`$, and $`F_\rho =F/\rho `$, $`F_\eta =F/\eta `$. The metric on the $`S^2`$ fiber of the twistor space of $`M`$ takes the form:

$$|d_A\stackrel{}{u}|^2=(du^i+ϵ^{ijk}A^ju^k)^2,$$
(13)
where $`\stackrel{}{u}=(u^1,u^2,u^3)`$ are constrained coordinates $`(_{i=1}^3(u^i)^2=1)`$ on the fiber and the connection $`A`$ is given by :

$$A^1=\frac{F_\eta }{F}d\rho +\left(\frac{1}{2\rho }+\frac{F_\rho }{F}\right)d\eta ,A^2=\frac{\sqrt{\rho }}{F}d\varphi ,A^3=\frac{\eta }{F\sqrt{\rho }}d\varphi +\frac{1}{F\sqrt{\rho }}d\psi .$$
(14)
Substituting (11) and (13) into (2) we obtain the $`G_2`$ metric of our models. As expected, this leads to rather complicated formulae, which can be found in  and will not be reproduced here.

As explained in , the singular loci of $`Y`$ project onto the boundary of $`^2`$. The latter is obtained by adding the point at infinity to the upper half plane: this results in a disk whose boundary corresponds to $`\rho =0`$. In agreement with the topological analysis of , the vertical loci project to points $`P_1\mathrm{}P_n`$ lying on this boundary, while the horizontal loci project to the segments $`P_kP_{k+1}`$ connecting them. The points $`P_k`$ correspond to $`\rho =0`$ and $`\eta =\eta _k`$, where $`\eta _k=\frac{\nu _k^1}{\nu _k^2}`$. The circular polygon (with vertices $`P_k`$) obtained in this manner can be identified topologically with the polygon $`\mathrm{\Delta }_M`$ defined as the quotient of the characteristic polygon $`\mathrm{\Delta }`$ (see (8)) through the point reflection $`\iota :aa`$ of the plane:

$$\mathrm{\Delta }_M=\mathrm{\Delta }/\iota .$$
(15)
The polygon $`\mathrm{\Delta }_M`$ has $`n`$ vertices. Opposite edges $`e`$ and $`e`$ of $`\mathrm{\Delta }`$ cover the same edge of $`\mathrm{\Delta }_M`$ (figure 2).

Starting with the purely geometric background $`ds_{11}^2=ds^2(𝐑^{3,1})+ds_{G_2}^2`$ of 11d supergravity (see (2)), one can reduce along one of the $`T^2`$ circles to obtain a IIA solution with non-constant dilaton and non-vanishing RR 1-form. Modulo a common modular transformation of the vectors $`\nu _j`$, it suffices to consider the reduction along the $`\varphi `$-circle. The explicit form of the resulting IIA solution can be found in  and will not be reproduced here.

Parameterizing the $`S^2`$ fibers of $`Y`$ by spherical coordinates:

$$u^3=\mathrm{cos}\theta ,u^1=\mathrm{sin}\theta \mathrm{cos}\chi ,u^2=\mathrm{sin}\theta \mathrm{sin}\chi $$
(16)
with $`\theta [0,\pi ],\chi [0,2\pi ]`$, one obtains coordinates $`(r,\rho ,\eta ,\chi ,\theta ,\psi )`$ for the internal 6-dimensional part of the resulting IIA background. For every fixed value of $`r`$, one has a compact 5-manifold $`𝒩^5`$ with coordinates $`(\rho ,\eta ,\chi ,\theta ,\psi )`$. This is the $`U(1)_\varphi `$ reduction of the twistor space $`Y`$. Under this reduction, the distinguished locus of the twistor space corresponds to the locus $`𝒩_D`$ in $`𝒩^5`$ defined by $`\rho =0,\theta =0,\pi `$ or $`\rho =0,\eta =\eta _k`$ (for some $`k`$). In the limit $`\rho 0`$, the angular part of the metric turns out to depend only on $`\theta `$ and $`\xi \nu _ϵ^1\chi +\psi `$, where $`\nu _ϵ^1`$ is the first component of a two-component vector $`\nu _ϵ`$ related to the edge of $`\mathrm{\Delta }_M`$ covered by the edges $`e,e`$ of $`\mathrm{\Delta }`$ ($`\nu _{\pm e}`$ of (9) equals $`\pm \nu _ϵ`$). Hence the locus $`\rho =0`$ in $`𝒩^5`$ is only three-dimensional, while the sublocus $`𝒩_D`$ has dimension two. The horizontal part $`\rho =0,\theta =0,\pi `$ of $`𝒩_D`$ turns out to be a union of $`2n`$ components $`𝒩_e`$ indexed by the edges of $`\mathrm{\Delta }`$. The vertical part consists of the components $`𝒩_k`$ defined by $`\rho =0`$ and $`\eta =\eta _k`$.

The $`\rho 0`$ asymptotics of the IIA fields (dilaton, RR 1-form and metric) was analyzed in detail in . Fixing a component $`𝒩_e`$ or $`𝒩_j`$, the IIA background turns out to be strongly coupled in its vicinity unless $`\nu _e^1=0`$ (resp. $`\nu _j^1=0`$), i.e. unless the direction $`\varphi `$ used for reduction corresponds to the direction $`\nu _e`$ (resp. $`\nu _j`$) in the Lie algebra $`𝐑^2`$ of the $`T^2`$ isometry group. Remember from the previous section that $`𝒩_e`$ (resp. $`𝒩_j`$) corresponds to enhanced gauge symmetry if $`\nu _e`$ (resp. $`\nu _j`$) is not primitive. To insure weak coupling along all such loci, one must be able to choose $`U(1)_\varphi `$ in the direction of all non-primitive $`\nu _e`$ and $`\nu _j`$. This requires that all these vectors be proportional (over the reals), which is the condition for a ‘good isometry’ mentioned above.

For simplicity, we shall restrict to the case when $`U(1)_\varphi `$ is good. In this situation, all interesting components of the locus $`𝒩_D`$ are weakly coupled. Moreover, the computations of  show that the interesting loci $`𝒩_e`$, $`𝒩_j`$ carry $`|\nu _e^2|`$ resp. $`|\nu _j^2|`$ units of RR flux. Since $`\nu _e^1=\nu _j^1=0`$, one has $`|\nu _e^2|=m_e`$ and $`|\nu _j^2|=m_j`$, where $`m_e=gcd(\nu _e^1,\nu _e^2)`$ and $`m_j=gcd(\nu _j^1,\nu _j^2)`$. This confirms the prediction of  reviewed in the previous section, and shows that the appearance of enhanced $`SU(m_e)`$ or $`SU(m_j)`$ symmetry at these loci is due to the presence of coincident D6-branes777The equivalence of the two mechanisms of enhancement of gauge symmetry - coinciding D-branes in type IIA string theory and M2-branes wrapping collapsing two-cycles at an $`ADE`$ singularity - was first pointed out in .. The geometry at $`\rho =0`$ is depicted in figure 3 for a model with $`n=3`$ which admits a good isometry (for example the $`W_{p,p,q}`$ models of ).

T-dualizing along the remaining $`U(1)`$ isometry (parameterized by $`\psi `$) one finds that the $`\chi `$-circle is always collapsed to a point in the limit $`\rho 0`$ . Analyzing the asymptotic behavior of the IIB fields enables the computation of the relevant NS-NS/R-R fluxes. When the $`U(1)`$ isometry along $`\varphi `$ is good, one finds that the corresponding IIB loci support delocalized five-branes, once again confirming the predictions of .

###### Acknowledgments.

The present work was supported by the Research Foundation under NSF grant PHY-0098527.


text:
# Potential of Mean Force between a Spherical Particle Suspended in a Nematic Liquid Crystal and a Substrate

## I Introduction

The alignment of nematic liquid crystals can be controlled by surface chemistry or topography on the scale of a few nanometers. When the uniform alignment is disturbed, this event can be detected through the appearance of multidomains in an optical microscope between cross-polarizers. One possible event that can trigger this transition is the binding of a protein or a virus to the surface, thereby providing the basis for development of sensors gupta98 ; vannelson02 . Recent experiments have in fact shown that such liquid crystal sensors are remarkably sensitive to the presence of minute amounts of protein on a substrate gupta98 ; vannelson02 .

Several experiments using colloids in bulk or confined liquid crystals have studied the influence of local disturbances in a nematic host poulin97 ; poulin98 ; gu00 . Colloidal particles in liquid crystals have also been studied by computer simulations billeter00 ; andrienko01 . One of the questions that arises from these experiments pertains to the nature of the liquid crystal mediated interaction between the colloids, or a colloidal particle and a surface.

This interaction can be quantified through a potential of mean force; molecular simulations provide a means to calculate it directly by fixing or constraining a reaction coordinate of interest. In the particular case of liquid crystals, however, conventional simulation techniques such as Monte Carlo (MC) or molecular dynamics (MD) in regular ensembles are limited in their ability to sample phase space. The system can be trapped in local energy minima and, for MD simulations, the time scales associated with the motion of the colloid can be prohibitively long.

The majority of available methods for measuring the potential of mean force either require prior knowledge of the free energy of a system roux94 ; engkvist96 ; vondele00 , or constrain the system to a limited range of the reaction coordinate straatsma92 ; sprik98 . A novel method is proposed in this work to calculate the potential of mean force which does not require prior knowledge of the free energy and permits efficient sampling. The proposed technique consists of an expanded ensemble in the reaction coordinate, implemented within the formalism of a recently proposed density-of-states sampling technique wang01a . Related methods have been proposed by de Oliveira et al. deoliveira96  and by Engkvist et al. engkvist96 . This method is applied here to determine the potential of mean force between a spherical colloidal particle suspended in a confined nematic liquid crystal and a wall.

## II Simulation Methods

### II.1 Potential of Mean Force

The potential of mean force (PMF) provides a measure of the effective difference in free energy between two states as a function of one or several “interesting” degrees of freedom. These degrees of freedom can be real coordinates of the system or a combination of them. They are often called reaction coordinate(s) $`\xi `$. An integration of the derivative of the Hamiltonian with respect to the coordinate of interest is performed along the reaction coordinate to provide a local estimate of the PMF, $`w`$, between two points $`\xi _1`$ and $`\xi _2`$:

$$w(\xi _2)w(\xi _1)=_{\xi _1}^{\xi _2}𝑑\xi ^{}\frac{H}{\xi }$$
(1)
The reaction coordinate $`\xi `$ need not be directly observable, in the case of chemical reactions, for example, it can be the extent of reaction, and the potential of mean force provides a measure of the activation energy. It can also be a physical coordinate, such as the distance between two solute molecules in a solvent—the PMF then reflects the work required to bring two solute particles from an infinite to a finite separation. Both examples involve significant energy barriers which would severely limit the effectiveness of conventional MD or MC simulations for its calculation.

The PMF is connected to $`p(\xi )`$ (the probability of finding the system of interest in a state corresponding to a particular value of $`\xi `$) through the reversible work theorem by chandler87

$$w(\xi )=kT\mathrm{ln}p(\xi ).$$
(2)
Since the probability of visiting high energy states is low, their sampling and consequently the estimate of an energy barrier is poor. In order to force a system to sample the desired portion of phase space, it is often advantageous to fix the reaction coordinate or to apply a biasing potential that changes the free energy in a way that samples the barrier more efficiently. The former approach has the shortcoming of requiring a large number of simulations, so as to obtain a reasonable resolution for subsequent thermodynamic integration. We concentrate on the latter approach, as it offers the possibility of considerable improvements in sampling efficiency and accuracy.

### II.2 Expanded Density of States Method (EDOS)

The idea of applying a biasing potential appears in many forms and names in the literature. Umbrella sampling torrie74 ; beutler94  or multicanonical sampling are examples of such a strategy berg92a . In these techniques, the potential of the system is artificially altered in a way to lower free energy barriers and permit more uniform sampling of the altered phase space. Histograms of thermodynamic quantities are subsequently reweighted to recover results corresponding to the original potential ferrenberg88 . Another way of traversing phase space more efficiently is through the use of expanded ensembles lyubartsev92  in which intermediate states are introduced. These intermediate states do not necessarily represent states of the physical system; their role is to facilitate transitions between states separated by large energy barriers lyubartsev92 ; escobedo95b . Expanded ensemble formalisms have the disadvantage that visits to a specific state are dictated through unknown weights which must be determined iteratively.

In recent years, new methods have been proposed to calculate the density of states of a system directly engkvist96 ; deoliveira96 ; wang01a . Here we refer to these techniques as density of states Monte Carlo methods (DSMC). These methods introduce a biasing potential, which is a running best estimate of the inverse of the density of states as a function of energy, $`g(E)`$. It is obtained by performing a random walk in energy space. The main attractive feature of these methods is that $`g(E)`$ need not be known a priori; it converges self-consistently to the correct density of states. DSMC has been successfully applied to such diverse systems as Potts and Ising models wang01b , Lennard-Jones systems yan02 , binary Lennard-Jones glasses faller01sc , proteins rathore02 , polymers jain02 , and bulk liquid crystals kim02sb .

In this work, the DSMC formalism is combined with an expanded ensemble where the expanded states correspond to a reaction coordinate. In order to span the entire range of expanded states and visit them with uniform probability during the simulation, the states are weighted accordingly. In a sense, the weighting factors act as a biasing potential. An idea similar in spirit was proposed earlier by Engkvist and Karlström engkvist96 ; their method, however, requires a specific adaptation of the way in which the weights are updated for any new system to be investigated. Our new technique does not suffer from this problem and therefore can be viewed as self-consistent.

Consider a system of $`N`$ particles in volume $`V`$ and at temperature $`T`$. The system is expanded in $`M`$ subsystems along a reaction coordinate $`\xi `$ such that the $`m^{th}`$ state corresponds to a particular value of $`\xi =\xi _m`$. For concreteness, the reaction coordinate is set to be a distance $`X`$ between two solute particles in a solvent. A set of $`M`$ states is defined by assigning to each of them a value of $`X`$ in the range $`[X_{m1},X_m]`$, with $`X_0`$ and $`X_{m=M}`$ being the absolute bounds of the region of interest. The partition function $`\mathrm{\Omega }`$ of the expanded ensemble is given by

$$\mathrm{\Omega }=\underset{m=1}{\overset{M}{}}Q(N,V,T,m)g_m=\underset{m=1}{\overset{M}{}}Q_mg_m,$$
(3)
where $`Q_m`$ and $`g_m`$ denote a canonical partition function and a positive weighting factor, respectively, each corresponding to a particular state $`m`$. During the course of the simulation solvent and solute particles are rearranged through arbitrary Monte Carlo moves, e.g., trial random displacements. Whenever one of the solute particles is displaced, the distance $`X`$ between two solutes assumes a certain value that corresponds to some state $`m`$. The probability $`p_m`$ with which a state $`m`$ is visited is related to the partition function through

$$p_m=\frac{g_mQ_m}{\mathrm{\Omega }}.$$
(4)
In general, for any two states $`m`$ and $`k`$, the following relationship holds:

$$\frac{Q_m}{Q_k}=\frac{p_mg_k}{p_kg_m}.$$
(5)
The free energy difference between states $`m`$ and $`k`$ can be inferred from the weighting factors and the population of the respective states according to

$$\beta [w(\xi _m)w(\xi _k)]=\mathrm{ln}\frac{Q_m}{Q_k}=[\mathrm{ln}g_k\mathrm{ln}g_m]+[\mathrm{ln}p_k\mathrm{ln}p_m].$$
(6)
The relative population of the states is only used as an indicator of uniform sampling. The physical meaning of the weighting factors is apparent in Eq. (4): if the states are sampled uniformly, the free energy difference between states $`m`$ and $`k`$ is given by the natural logarithm of the ratio of the corresponding weighting factors. This stresses the importance of finding an efficient scheme that allows the weights to converge to the their “true” free energy value rapidly and accurately.

Continuing with the example of two solute particles in a solvent, a transition from an old to a new state occurs when the distance $`X`$ between the two solutes changes accordingly (molecular rearrangement of the solvent particles does not change the system’s state). The criteria for accepting a trial move from an old to a new state are given by

$$P_{\text{acc}}(\mathrm{𝑜𝑙𝑑}\mathrm{𝑛𝑒𝑤})=\mathrm{min}\{1,\mathrm{exp}[\beta (U_{\text{new}}U_{\text{old}})(\mathrm{ln}g_{\text{new}}\mathrm{ln}g_{\text{old}})]\}=\mathrm{min}\{1,\mathrm{exp}[\beta \mathrm{\Delta }U\mathrm{\Delta }\mathrm{ln}g]\}.$$
(7)
Every time a state is visited, the corresponding weight is modified by a convergence factor, $`f`$:

$$gfg$$
(8)
and the histogram of that state is updated. Note that in practice it is more convenient to work with $`\mathrm{ln}g`$, rather than $`g`$ itself since the former is used directly in the acceptance criteria; hence $`\mathrm{ln}f`$ is added to $`\mathrm{ln}g`$. According to Eqs. (7) and (8), if the original state has a higher energy than the new one, the new state will be visited more easily because $`\mathrm{exp}(\beta [U_{new}U_{old}])`$ is larger than unity and $`g_{\text{new}}`$ will be updated more frequently than $`g_{\text{old}}`$. Eventually, $`g_{\text{new}}/g_{\text{old}}`$ becomes smaller than unity and forces visits to the old, higher energy (and less visited) state. By construction, if $`g`$ is underestimated, i.e. if the contribution of the state to the total partition function is too small, the system can readily visit the under-visited state, and its weight is increased accordingly until a balance is reached; the converse also holds. Once a global histogram is sufficiently flat (the minimum population is enforced to be at least 85% of the average), a cycle is assumed to be “complete” and the convergence factor is altered according to an arbitrary, monotonically decreasing function; in the present implementation we use the square root. The random walk cycle is then resumed. The flatness condition stems from the fact that, if the inverse density of states were used as a weighting function, the histogram of energy would be perfectly flat, as the probability of visiting a state would be uniform.

Several issues are worth pointing out. First, in previous applications of DSMC, simulations proceeded until $`\mathrm{ln}f`$ reached a low threshold value; $`\mathrm{ln}f<10^8`$ has been commonly used. Such precision was important as $`g`$ reflected the density of states of the entire system over many orders of magnitude in energy. In the application discussed here, such precision is not required; the range of free energy or $`\mathrm{ln}g`$ is much narrower (temperature is fixed), and only a single or several reaction coordinates are of interest. Second, as can be inferred from the description of the method, the weights involved in the acceptance criteria are updated continuously, and the method does not obey detailed balance. Once the weights have been determined, however, the simulation can be rerun in a regular expanded ensemble without updating the weights, to verify that state population is indeed uniform. Third, the weight updating scheme proposed here may be used as an efficient tool for adequate sampling of any system with inherent, high energy barriers. If one is interested in a particular quantity that can be measured directly during the simulation (e.g. the mean force), it is sufficient to stop the simulation once good sampling of that quantity is achieved. Here we are interested in the potential of mean force; once a good estimate is obtained, the mean force can be integrated to yield the potential of mean force.

### II.3 Umbrella Sampling

In order to compare the results of the proposed method to more traditional techniques, umbrella sampling simulations were conducted over a sub-range of the reaction coordinate of interest. The particular implementation of the umbrella potential method employed here is discussed in what follows. As before, an artificial potential is introduced to prevent the system from being trapped in local energy minima. If multiple energy minima are present, i.e. the system is characterized by a rough energy profile, it is more effective to subdivide the entire range of the reaction coordinate into smaller overlapping ranges or “windows”, and to perform multiple simulations for each window. If the reaction coordinate is a distance, as in our case, a natural choice for the potential is a harmonic spring potential. The rapid increase at some distance from the spring equilibrium position guarantees that the system is contained within a chosen distance window. Alternatively, instead of applying a pre-determined umbrella potential, one can iteratively arrive at an umbrella potential that should be close to $`w(\xi )`$ through a multicanonical simulation. To this end one would record the histogram using a given umbrella potential, Boltzmann invert it to get an estimate of the relative free energy, and add this estimate to the potential faller02c .

For any given window, the PMF is defined as

$$\beta w(\xi )=\mathrm{ln}\left\{\frac{\mathrm{}e^{\beta U}𝑑r_{n+1}\mathrm{}𝑑r_N}{\mathrm{}e^{\beta U}𝑑r_1\mathrm{}𝑑r_N}\right\}=\mathrm{ln}\left\{\frac{Z(\xi )}{Z}\right\}$$
(9)
where $`\xi `$ is a $`n`$-dimensional set of reaction coordinates related to coordinates $`r_1\mathrm{}r_n`$. If the umbrella potential $`U_{\text{umb}}`$ is introduced then the PMF can be computed from:

$$\beta w(\xi )=\mathrm{ln}Z_{\text{umb}}(\xi )\beta U_{\text{umb}}(\xi )\mathrm{ln}Z$$
(10)
Here $`Z_{\text{umb}}`$ and $`Z`$ correspond to the partition functions of the system with and without an umbrella potential, respectively. Another way to express the PMF is

$$\beta w(\xi )=\mathrm{ln}g_{\text{umb}}(\xi )\beta U_{\text{umb}}\mathrm{ln}\frac{Z_{\text{umb}}}{Z},$$
(11)
where $`g_{\text{umb}}(\xi )`$ denotes the correlation function or the probability of finding a system at a particular value of $`\xi `$. Note that the last term in the above equation is a constant (which is different for every window). Since we are interested in the relative free energy and it is a continuous function of $`\xi `$, the overlapping ranges of the curves $`w(\xi )`$ have to be collapsed onto one master curve by fitting an offset for each window. Once again, the disadvantages of the method in this rendering include uncertainty in identifying an optimal umbrella potential, in choosing the number of necessary windows, and in fitting the final results for complex systems.

## III Simulation Details

The system considered in this work comprises 11460 liquid crystal particles confined between two soft repulsive walls at $`z=0`$ and $`z=Z_{\text{wall}}`$ (see Figure 1a). Simulations were conducted at a constant temperature $`T^{}=1.0`$ and an average number density of $`\rho ^{}=0.335`$. All variables are given in reduced units, e.g. $`\rho ^{}=N\sigma _0^3/V`$, $`T^{}=kT/ϵ_0`$, where $`\sigma _0`$ and $`ϵ_0`$ are length and energy parameters which are set to unity. The shape of the simulation box is orthorhombic; the side lengths are equal in the $`x`$ and $`y`$ directions. Periodic boundary conditions apply in $`x`$ and $`y`$. Mesogens $`i`$ and $`j`$ interact via a shifted and truncated Gay–Berne (GB) potential gay81 :

$`U_{ij}`$ $`=`$ $`\{\begin{array}{cc}4ϵ_0(\varrho _{ij}^{12}\varrho _{ij}^6)+ϵ_0,\hfill & \varrho _{ij}^6<2\hfill \\ 0,\hfill & \varrho _{ij}^6>2\hfill \end{array}`$ (12)
$`\varrho _{ij}`$ $`=`$ $`(|𝒓_{ij}|\sigma _{ij}+\sigma _0)/\sigma _0`$ (13)
$`\sigma _{ij}`$ $`=`$ $`\sigma _o\left[1{\displaystyle \frac{\chi }{2}}\left\{{\displaystyle \frac{\left(\widehat{𝒓}_{ij}\widehat{𝒖}_i+\widehat{𝒓}_{ij}\widehat{𝒖}_j\right)^2}{1+\chi \widehat{𝒖}_i\widehat{𝒖}_j}}+{\displaystyle \frac{\left(\widehat{𝒓}_{ij}\widehat{𝒖}_i+\widehat{𝒓}_{ij}\widehat{𝒖}_j\right)^2}{1\chi \widehat{𝒖}_i\widehat{𝒖}_j}}\right\}\right]^{1/2}`$ (14)
$`\chi `$ $`=`$ $`{\displaystyle \frac{\kappa ^21}{\kappa ^2+1}}.`$ (15)
where $`\widehat{𝒖}_i`$ and $`𝒓_{ij}`$ denote molecular orientations (along the main axis of the ellipsoid) and intermolecular vectors ($`𝒓_{ij}=𝒓_i𝒓_j`$, where $`𝒓_i`$ is the position of particle $`i`$); unit vectors are identified by hats ($`\widehat{𝒓}_{ij}=𝒓_{ij}/|𝒓_{ij}|`$). For details see Figure 1b. Parameter $`\kappa =3`$ in this work is the length to width ($`\sigma _0`$) ratio. The interactions of a GB mesogen $`i`$ with all interfaces are described by the same potential (12), where the $`\varrho _{ij}`$ are defined as follows for a sphere and a surface:

$`\varrho _{i\text{sph}}`$ $`=`$ $`\left(|𝒓_i𝒓_{\text{sph}}|R+\sigma _0/2\right)/\sigma _0,`$ (16)
$`\varrho _{i\text{wall}}`$ $`=`$ $`\left(|𝒛_i|+\sigma _o/2\right)/\sigma _0.`$ (17)
In this case, $`𝒓_{\text{sph}}`$, $`𝒓_i`$, and $`|𝒛_i|`$ denote the position vectors for the centers of mass of the sphere, mesogen $`i`$, and the normal distance from a surface, respectively; $`R`$ represents the radius of the sphere. This potential has been employed in a previous study of a colloidal particle in a bulk liquid crystal andrienko01 , and is known to reproduce topological defects around the particle that are in agreement with theory and experiment. The force on a sphere is calculated as the derivative of the interaction potential between the liquid crystal and the sphere with respect to $`z`$; due to symmetry in the $`(x,y)`$-plane, only the $`z`$-component of the force contributes. The surface and sphere interact as hard bodies, i.e. the minimum distance possible between a surface and the center of mass of the sphere is the radius of the sphere, which is set to 3.

Before proceeding with this system, we characterized the confined Gay-Berne system without a colloid in independent simulations. The separation between the walls was chosen to be $`Z_{\text{wall}}=34`$. A second rank order parameter $`P_2`$ is defined as

$$P_2=\frac{1}{2}3(\widehat{𝒖}_i\widehat{𝑷})^21,$$
(18)
where $`\widehat{𝑷}`$ is the average orientation of the system. Figure 2 shows the resulting profiles of density and order parameter. One can see that several layers are formed at the surface. In the midsection between the walls both density and $`P_2`$ profiles are almost flat. In that midsection, a sphere is unlikely to “feel” the presence of a substrate. On the other hand, the movement of the sphere in the vicinity of the surface may be severely limited or dictated by the layered structure of the liquid crystals.

## IV Results and Discussion

The quantity of interest is the free energy profile for a colloidal sphere in a confined nematic liquid crystal as a function of distance between the center of mass of the sphere and the surface. From the results shown in Figure 2 the fluid structure is rich near the surface and bulk-like for $`z>13`$. Hence, we concentrate on the region $`3z_{\text{sph}}13`$, where $`z_{\text{sph}}`$ denotes the position of the center of mass of a sphere and plays the role of a reaction coordinate. The lower bound is the radius of the sphere. As a first trial, $`NVT`$ MC simulations were performed with $`z_{\text{sph}}`$ fixed at specific values. A plot

(Figure 3) of the mean force acting on a sphere from these calculations indicates that indeed the effective force from the wall on a sphere is increasingly repulsive (positive) up to $`z_{\text{sph}}=3.8`$ and then decreases to zero around the value of 4. The force then changes from repulsive to attractive and back to repulsive and becomes seemingly constant and neutral for $`z_{\text{sph}}>9`$. From the graph we can infer that at $`z_{\text{sph}}=4`$ the system reaches a stable equilibrium. This is confirmed by simulations of a freely moving sphere; from its original position of $`z_{\text{sph}}=3`$, the sphere travels within a few MC cycles to a position of $`z_{\text{sph}}=4`$, where it undergoes only slight oscillations. The structure of the liquid crystals near the surface influences considerably the mean force on a sphere, but conventional $`NVT`$ Monte Carlo simulations are unable to elucidate the details of these interactions.

Next, we performed umbrella sampling simulations for $`z_{\text{sph}}`$ near the surface. In that region, the energy profile was found to vary strongly. Highly stiff springs are necessary, which in turn result in very narrow windows. To widen the range of the windows and improve their overlap, we applied spring potentials of the form $`U_{\text{umb}}=A(z_{\text{sph}}z_0)^4+B(z_{\text{sph}}z_0)^2`$, where $`z_0`$ is an equilibrium position for a spring and $`A>0`$ and $`B<0`$. Such potentials have a “W” shape and therefore result in a slight dispersion of a sphere close to $`z_0`$. Depending on the behavior of the system, at some points we had to use hybrid springs that combined a stiff harmonic spring on one side of $`z_0`$ and a dispersing potential on the other. As a result, 13 windows were necessary to explore the narrow range of $`z`$ between $`[3.0,4.2]`$. The corresponding umbrella potentials are given in Table 1.

Figure 4 shows corresponding histograms collected over a simulation of $`3`$$`5\times 10^5`$ MC cycles from the initially equilibrated configurations. Not surprisingly, a histogram centered at 4 is wide and Gaussian in appearance, whereas the others are sharp and narrow. Histogram data $`g_i(z_{\text{sph}})`$ are transformed by taking the natural logarithm and subtracting the corresponding umbrella potentials. Even though each data set is offset by a constant, there appear to be at least two separate regions where different functions must be fitted. We fit the first 12 sets and the last set separately using third degree polynomials by a least squares procedure modified to include constant offsets (see Figure 5); the third degree was chosen as the minimum required to have an inflection point, i.e. the derivative of the force is zero. The intersection of forces (negative derivatives) is taken as the matching point; the curves of the potential of mean force are joined at that point $`(z_{\text{sph}}=3.8,F(z_{\text{sph}}=3.8)=164)`$ setting the PMF at $`z_{\text{sph}}=3`$ to zero. The resulting force profile (see Fig. 6) is in agreement with results of canonical simulations but reveals more detail.

Next, the application of the new method is demonstrated on the same system. A state $`m`$ is defined by the position of the center of mass of a sphere in the range $`(m1)\mathrm{\Delta }z<z_{\text{sph},m}m\mathrm{\Delta }z`$. The entire range of $`[3,13]`$ is divided into nine overlapping windows of width 2: the first window covers $`3z_{\text{sph}}<5`$, the second window covers $`4z_{\text{sph}}<6`$, etc. Each window consists of 200 states with $`\mathrm{\Delta }z=0.01`$. A move to a new state, i.e. a move of a sphere in the $`z`$ direction (its $`x`$ and $`y`$ coordinates are kept at zero throughout the entire simulation) is attempted every Monte Carlo cycle; a cycle consists of a trial displacement or rotation of all liquid crystal particles. Initially, all weighting factors $`g_m`$ and the convergence factor $`f`$ are given a value of $`e^{0.1}`$. Simulations proceed for three DSMC cycles: that is, $`\mathrm{ln}f`$ is halved twice. The resulting mean force and relative weights are averaged between windows in the area of overlap. Figure 6b shows the results for the mean force. The corresponding free energy profiles are shown in Figure 7 (integrated mean force; PMF was arbitrarily set to 0 at $`z_{\text{sph}}=3)`$. It is evident from the above plots that the EDOS method provides improved sampling compared to simple $`NVT`$ or umbrella potential simulations. Moreover, more details of the interaction of the colloidal sphere with the liquid crystals and the surface are visible (compare Figures 3 and 6a). Note that with only 9 windows it is possible to cover a substantially wider range than with 13 windows in umbrella sampling.

It is apparent that the mean force on the sphere has different characteristics in various regions in proximity of the surface. Close to the surface, the force is strongly repulsive and reaches its maximum ($`F_{max}170`$) at $`z_{\text{sph}}3.8`$. The force then decreases rapidly, dropping to zero at $`z_{\text{sph}}=4`$ and eventually reaching a minimum $`(F_{min}35)`$ at $`z_{\text{sph}}=4.3`$. In the interval $`4<z_{\text{sph}}<11`$ the force exhibits a damped oscillatory behavior with a wavelength of about 2.5–2.6. The force oscillations die out within about 2–3 of these wavelengths. At that point ($`z_{\text{sph}}11`$) the force becomes effectively zero.

Clearly, the layered structure of the liquid crystals has a strong effect on the behavior of a suspended colloidal sphere. The effect appears to be most pronounced in the immediate vicinity of the surface, where the density and the order parameter are comparable to those of a smectic or even solid layer. Figure 8 shows pictures of the colloid seen from below (the surface was removed for clarity). The sphere severely disrupts the first, dense layer of liquid crystals; as the sphere moves away from the surface, the cavity caused by it gradually fills up. At $`z_{\text{sph}}=3.8`$ the layer is almost completely restored. At $`z_{\text{sph}}=4`$, the first layer is fully repaired; this point corresponds to a stable state of the system as the force on the sphere vanishes (see above).

As for subsequent undulations in the force profile, their wavelength approximately coincides with the thickness of the liquid crystalline layers close to the surface. If we superimpose the density profile of the confined system without the sphere (see Figure 9) onto the PMF curve, the oscillations are almost in phase. The PMF curve is slightly shifted to larger distances from the surface (offset in $`z_{\text{sph}}0.25)`$. The minima in the PMF practically match those of the density profile; our results indicate that the sphere prefers to reside between adjacent LC layers.

## V Conclusion

A new method is proposed in which a density of states Monte Carlo formalism is used to self-consistently calculate the weights for an expanded ensemble technique. This method permits efficient calculation of the potential of mean force of a complex system. Its application is illustrated in the context of a colloidal sphere dispersed in a confined nematic host.

The method allows one to determine with high resolution the mean force on the colloid as a function of its separation from the surface; the results agree qualitatively with less efficient methods such as umbrella sampling and conventional canonical Monte Carlo. Clearly, the layered structure of the liquid crystals near the surface profoundly influences the mean force on the colloid. The snapshots of the system at several distances indicate that the disruption of the first surface layer is most likely responsible for the strong repulsive force very close to the surface. Once the “hole” caused by the sphere is completely filled in, the force decreases to neutral at $`z_{\text{sph}}=4`$, resulting in a preferred position for the colloid reflected by a minimum in its free energy (potential of mean force). At larger distances, both the force and the potential of mean force profiles display undulations that are almost in phase with the oscillations in the density profile of the confined system. It is worth pointing out again that all intermolecular interactions in the simulation are strictly repulsive. The attractive regions in the potential of mean force between the colloid and the surface have therefore purely entropic origins. Additional studies are currently under way to elucidate the nature of the interaction between the surface and the colloid for different sphere sizes and different ancoring conditions.

## Acknowledgments

This work was supported by the University of Wisconsin MRSEC for nanocomposites. RF also thanks the DFG (German Research Foundation) for financial support through the Emmy–Noether program.