format: text

text:
# Decoherence of electron spin qubits in Si-based quantum computers

## I Introduction

The prospect of quantum computing (QC) has caused great excitement in condensed matter physics. If a set of qubits can be maintained in a coherent, controllable many-body state, certain very difficult computational problems become tractable. In particular, successful QC would mean a revolution in the areas of cryptography shor  and data-base searching grover . In addition, it would mean a great advance in general technical capabilities, since the control of individual quantum systems and their interactions would represent a new era in nanotechnology.

However, from a practical point of view, a dilemma presents itself immediately. On the one hand, one wishes to control quantum degrees of freedom using external influences, since that is how a quantum algorithm is implemented, and to measure them, since that is the output step. On the other hand, the system must be isolated from the environment, since random perturbations will destroy the quantum coherence that is the whole advantage of QC. This is the isolation-control dilemma, and it leads to a very rough figure-of-merit $`F`$ for any quantum computer. If we define the decoherence time $`\tau `$ as the time it takes to lose quantum coherence, and the clock speed $`s`$ (roughly the inverse of the time to run a logic gate), then the figure of merit is $`F=s\tau `$ and a practical machine should satisfy $`F>10^410^5`$ at least. If the clock speed is limited only by standard electronics, then we may be able to achieve $`s10^9Hz`$. This would imply that $`\tau =1ms`$ is a lower limit for the decoherence time.

The dilemma has not yet been solved, though a number of solutions have been proposed. A particularly attractive solution is to use spin degrees of freedom as qubits. Nuclear spins interact relatively weakly with their environment because the coupling, proportional to the magnetic moment, is small. Yet there has grown up a sophisticated technology (NMR) for the manipulation of the spins, and some rudimentary computations have been performed chuang . Readout is the main difficulty with this approach, since the field created by a single moment is tiny, and pure states cannot be achieved in the macroscopic samples used.

Electron spins also interact weakly with the environment in some circumstances: relaxation times in excess of $`10^3s`$ have been measured for donor bound states of phosphorus-doped silicon (Si:P) feher . The corresponding manipulation technology (ESR) has also reached a high level of sophistication, but the magnetic moment of the electron exceeds that of the nucleus by 3 orders of magnitude, which presents problems of isolation. Readout should be easier than for nuclei, since detection of single electron charges is certainly possible kastner , and spin-dependent detection is not far out of the reach of current technology.

Solid-state implementations of QC are particularly attractive because of the possibility of using existing computer technology to scale small numbers of qubits up to the $`10^5`$ or so that would be needed for nontrivial computations. The first paper to propose using the electron spin in a quantum dot subjected to a strong DC magnetic field was that of Loss and DiVincenzo loss . Kane kane  proposed employing the nuclear spin in the Si:P system as the qubit. A specific structure consisting of silicon-germanium (SiGe) layers was proposed by Vrijen et al. vrijen . This structure incorporates the idea that the $`g`$-factor of an electron can be changed by moving it in a Ge concentration gradient, allowing individual electron to be addressed by the external AC field. A different SiGe structure has been proposed by Friesen et al. eriksson . This structure is designed so that the electron number on the dots, and the coupling between the dots, can be carefully controlled.

Solid-state implementations must also face the isolation-control dilemma. Decoherence times must exceed the $`1ms`$ number in the actual physical structures that are needed for the operation of quantum algorithms. In this paper, we examine whether this can be the case for some of the existing proposals based on electron-spin qubits. In the process, we hope to learn something about modifications to these structures that can increase $`\tau `$. We shall focus on low temperature operation, since, as we shall see, this will probably be necessary in order to obtain sufficiently large $`\tau `$.

We can build on a large body of work, both theoretical and experimental, from the 1950’s and 1960’s on ESR in doped semiconductors. In a series of papers, Feher and collaborators feher  investigated the relaxation time for the spin of electrons bound on donor sites in lightly doped Si. At sufficiently low temperatures, the relaxation time $`T_1`$ is dominated by single-phonon emission and absorption. In the presence of spin-orbit coupling (SOC), this can relax the spin, causing decoherence. The theory was worked out by Hasegawa hasegawa  and Roth roth .

It must of course be recognized that this spin-lattice relaxation time is not necessarily to be identified with the decoherence time. The decoherence time is the shortest time for any process to permanently erase the phase information in the wavefunction. This may mean the phase for a single spin, but it also means that the relative phases of the wavefunctions of different spins must also be preserved, so that processes that cause mutual decoherence must also be taken into account. The actual decoherence time is the minimum of all of these times. A spin relaxation time in excess of $`1ms`$ is a necessary, not a sufficient, condition for the viability of a solid-state electron-spin QC proposal.

A QC must have precise input as well as an accurate algorithm. Preparation of the spin state is often proposed to be done by thermalization of the spin system at a low temperature. The time to do this actually sets an upper limit on the relaxation time of whatever processes thermalize the spins to the lattice. A limit of perhaps $`110`$ $`s`$ is a reasonable requirement.

This paper focuses on $`T_1`$, the time for relaxation of the longitudinal component of the magnetization, by spin-phonon interactions. These processes cause real spin-flip transitions. They occur at random times and thus indubitably cause decoherence. In addition to these processes characterized fully by $`T_1`$, there are processes which introduce random phase changes in the spin wavefunctions. To characterize all such processes by a single ”dephasing time” $`T_2`$ will usually not be sufficient for understanding the operation of a multi-qubit system hu1 . Difficulties of definition arise, and care must be taken to specify which phase is involved and to what extent it is randomized. For a single spin system there is no ambiguity. The $`2\times 2`$ density matrix $`\rho _{ij\text{ }}`$for the qubit with cylindrical symmetry involves only 2 independent parameters $`\rho _{11}\rho _{22}`$, and $`\rho _{12}.`$ The time dependence of $`\rho _{11}\rho _{22}`$ after a system preparation is exponential with decay constant $`T_1`$ and represents the return of the longitudinal component of the magnetization to its equilibrium value. The decay is due to inelastic transitions of the type calculated in this paper. $`\rho _{12}`$ on the other hand, is nonzero only if the preparation of the spin state has a transverse component: $`S_x(t=0)0.`$ The decay of this quantity represents the irreversible conversion of this state to an incoherent mixure of ”up” and ”down” states. Again, this is genuine decoherence of the spin state, since the phase information cannot be recovered. The time dependence of $`\rho _{12}`$ when the spin is in a strong field has been calculated by Mozyrsky et al. mozyrsky  using a Markovian approximate master equation. They find that the time scale of the decay due to spin-phonon coupling is very short, of the order of the time for a phonon to cross the electron’s wavefunction, which is about $`10^{10}s`$. But the decay is incomplete, with $`\rho _{12}`$ retaining all but $`10^8`$ of its original value. Their calculation was for Si:P, but a very similar result should hold for the dot case. Decoherence at this level is certainly acceptable for quantum computing. These authors also pointed out that the decay of the remainder of $`\rho _{12}`$ is due to the spin-flip processes computed in this paper. If $`T_2`$ is defined as the dominant decay time of the off-diagonal density matrix element, then $`T_1=T_2`$ for spin-phonon processes.

A quite different source of decoherence is the hyperfine coupling to nuclear spins. The nuclear spins produce an effective random magnetic field on the electrons. Very recent calculations using semiclassical averaging techniques efros  have obtained a very short relaxation time $`T_\mathrm{\Delta }1ns`$ for GaAs-based dot systems in a strong field. This represents the decay of the transverse magnetization of an ensemble of dots. This is a dephasing time, but not a decoherence time. The electrons spins precess in what is effectively the frozen field of the nuclei. This field is spatially random, and the differential precession of the electron spins leads to the magnetization decay. However, this is not an irreversible loss of the phase information of the collective wavefunction. Spin echo exeriments are a very beautiful demonstration of precisely this point. This ”inhomogeneous broadening” presents challenges for the calibration and operation of quantum computers, but does not destroy coherence.

Finally, in any implementation based on electron spin qubits, there will certainly exist small interactions between the spins themselves. The dipole-dipole interaction, for one, cannot be avoided, and there may be indirect spin-spin interactions mediated by the gates. A recent paper suggests that these interactions set the fundamental time scale $`T_M`$ for Si quantum dot implementations of QC desousa . These interactions do produce experimental broadening of ESR lines in experiments on bulk systems, and this might be taken as decoherence. In our view, however, these interactions do not destroy the coherence of a state. The system is the set of all the qubits. During the course of a quantum algorithm they are collectively in a pure state (in principle). Any decoherence that destroys the purity of the state comes from averaging over the unknown states of the environment. The broadening that comes from dipole-dipole interactions comes, in NMR and ESR calculations, from averaging over the states of the system itself, which is not an appropriate method for calculating decoherence. The effect of qubit-qubit interactions that cannot be turned off is to complicate the quantum algorithm. A quantum algorithm is a unitary transformation that must always include the effect of the system Hamiltonian (including dipole-dipole interactions) in addition to external operations. In every case except for very simple ones, this algorithm must first be computed, presumably with the help of a classical computer. This step in QC may be termed ”quantum compilation”. The issue that qubit-qubit interactions raise is not one of decoherence, but rather whether the determination of the algorithm, the complation step, becomes prohibitively difficult. This could happen for two reasons. One is that the interactions are so poorly known that they cannot be corrected for. It seems likely that quantum error correction can resolve this difficulty. A second and more interesting possibility that the interactions convert the computation of the algorithm itself into a problem that grows exponentially with the size of the system. We regard this as an open question and a deep one, that combines many-body theory with algorithm design and error correction. We note that in NMR implementations the interaction between the qubits also cannot be turned off, but it can be canceled by refocusing chuang2 .

In this paper, our aim is to evaluate the importance of spin-phonon coupling as a source of decoherence in quantum dot qubits. Fundamentally, the issue is whether the long relaxation times $`T_1`$ observed at low temperatures in bulk Si:P carry over to SiGe dots proposed for QC. In Sec. II we introduce the structures that we are interested in. In Sec. III the physics of the spin-phonon relaxation mechanism is described qualitatively. In Sec. IV calculations and results for pure Si systems are given as a function of the critical design and operating parameters. In Sec. V are found the correponding results for structures containing Ge. Sec. VI is the conclusion.

## II The strained silicon quantum well

Si-Ge heterostructures are utilized widely in the digital electronics industry, and presently have the shortest switching times of any device. One reason for their success lies in the ability to engineer structures of near perfect purity, with control over thicknesses and interfaces that approaches atomic precision – a technological tour-de-force. An equally key achievement has been the harnessing of strain as a tool to control band offsets in heterostructure devices. This paper presents calculations of spin relaxation for real SiGe structures such as those proposed by Vrijen et al. vrijen  and Friesen et al. eriksson . Accordingly, we have calculated the electron wavefunctions in quantum wells, which is needed as input for these calculations. Details of these calculations are presented in Ref. eriksson , and will not be repeated here. In this section we only describe those aspects of the calculations that are germane to spin relaxation.

Quantum wells are constructed by sandwiching a very thin layer of one material between two others. Electrons can be confined in the quantum well layer when the conduction band offsets produce a potential well. The key to this technology is therefore to understand the band structure of the various layers. In this section we will consider a particular class of wells formed of pure Si, sandwiched between barrier layers of SiGe. We will find that this is optimal from the standpoint of spin coherence. Metallic gates or impurities create zero-dimensional bound states that define a quantum dot. We first review briefly effective mass theory for dots in pure, unstrained Si, then unstrained SiGe, and finally strained Si.

In pure Si, the $`\mathrm{\Delta }`$ conduction band minima occur near the symmetry points $`X`$, in the directions $`\{001\}`$. In a perfect Si crystal these minima are six-fold degenerate, the valleys being equivalent. In the dot the electron feels a potential $`V_g(\stackrel{}{r})`$ in addition to the atomic potential, which lifts the degeneracy, though the splittings are not large. The spatial variation of $`V_g(\stackrel{}{r})`$ is on length scales generally much longer than the lattice spacing. For the moment, we shall assume that the electron is in the ground state of $`V_g(\stackrel{}{r})`$ and ignore mixing with any excited states. To the extent that the scale of variation of $`V_g(\stackrel{}{r})`$ is much longer than the atomic spacing there are six nearly degenerate ground states. This is referred to as the ”valley degeneracy”. The wavefunctions can be written askohn57

$$\mathrm{\Phi }_n(\stackrel{}{r})=\underset{j=1}{\overset{6}{}}\alpha _n^{(j)}F_j(\stackrel{}{r})\varphi _j(\stackrel{}{r}).$$
(1)
Here, $`\varphi _j`$ is a Bloch function of the form

$$\varphi _j(\stackrel{}{r})=u_j(\stackrel{}{r})e^{i\stackrel{}{k}_j\stackrel{}{r}},$$
(2)
where $`\stackrel{}{k}_j`$ are the six $`\mathrm{\Delta }`$ minima $`\{+k_0\widehat{x},k_0\widehat{x},+k_0\widehat{y},k_0\widehat{y},+k_0\widehat{z},k_0\widehat{z}\}`$ (we shall always use this ordering), and $`u_j(\stackrel{}{r})`$ are periodic functions with the same periodicity as the crystal potential $`V_p(\stackrel{}{r})`$. The $`F_{\pm z}`$ are envelope functions that satisfy the Schrödinger-like equation

$$[\frac{\mathrm{}^2}{2m_l}\frac{^2}{z^2}\frac{\mathrm{}^2}{2m_t}(\frac{^2}{x^2}+\frac{^2}{y^2})+V_g(\stackrel{}{r})]F_{\pm z}(\stackrel{}{r}))=(EE_{k_z}^{(\mathrm{\Delta })})F_{\pm z}(\stackrel{}{r}),$$
(3)
and are independently normalized to unity, similar to wavefunctions. Analogous equations can be given for the $`\pm \widehat{𝐱}`$ and $`\pm \widehat{𝐲}`$ minima. We see that $`F_x=F_x`$, $`F_y=F_y`$, and $`F_z=F_z`$, so only three independent envelope functions must be computed. $`m_l`$ and $`m_t`$ are the longitudinal and transverse effective masses associated with the anisotropic conduction band valleys. $`E_{k_z}^{(\mathrm{\Delta })}`$ is the $`\mathrm{\Delta }`$ conduction band edge at $`𝐤_z`$. The splitting of the degeneracy comes from corrections to this envelope fuction approximation. Different choices of the constants $`\alpha _n^{(j)}`$ determine the six states. Their values will be discussed in Sec. III. This formalism is a good approximation for both dot and impurity bound states, as the valley splitings are much smaller than the energy scales in 3.

Germanium is completely miscible in Si, forming a random alloy. For variable Ge content $`x`$, Si1-xGex exhibits materials properties that vary gradually over the composition range. The alloy lattice constant, $`a_0(x)`$, follows a linear interpolation between pure Si and Ge, known as Vegard’s law, quite accurately for all $`x`$: $`a_0(x)=(1x)a_{\mathrm{Si}}+xa_{\mathrm{Ge}}`$ windl98 . Electronic properties show an abrupt change in behavior near $`x0.85`$, where the Si-like $`\mathrm{\Delta }`$ minima cross over to four-fold degenerate, Ge-like, $`L`$ minima. In this work we focus on the range $`x0.5`$, which is strictly Si-like, though we will have some remarks below on Ge-rich structures. Throughout this range, properties such as effective mass and the dielectric constant vary only slightly from pure Si values. For our calculations, the most important parameter is the conduction band edge, $`E^{(\mathrm{\Delta })}(x)`$, which remains 6-fold degenerate in the range $`x0.5`$. The theory of the variation of $`E^{(\mathrm{\Delta })}`$with $`x`$ is not germane to the present work, and we simply quote the empirical result, linear in $`x`$, which is consistent with Ref. rieger93 :

$$\mathrm{\Delta }E^{(\mathrm{\Delta })}(x)=E^{(\mathrm{\Delta })}(x)E^{(\mathrm{\Delta })}(0)0.23x(\mathrm{eV}).$$
(4)
(We note, however, that Ref. schaeffler97  suggests a slope for $`E^{(\mathrm{\Delta })}(x)`$ of opposite sign.) The relatively weak variations of the effective mass and the dielectric constant will be ignored here.

We consider thin Si wells in which the Si layer grows pseudomorphically. The in-plane lattice constant, $`a_{}`$, must be the same for all the layers, causing a tetragonal distortion in the strained layer(s). Here we consider the case of strained Si grown on the (001) surface of relaxed Si1-xGex. The in-plane Si lattice constant depends on $`x`$ as

$$a_{}(x)=(1x)a_{\mathrm{Si}}+xa_{\mathrm{Ge}}.$$
(5)
Since $`a_{\mathrm{Ge}}>a_{\mathrm{Si}}`$, the Si is under tensile strain in the plane. Hence, the out-of-plane Si lattice constant, $`a_{}`$, is reduced according to continuum elastic theory:

$$a_{}(x)=a_{\mathrm{Si}}\left[12\frac{c_{12}}{c_{11}}\frac{a_{}(x)a_{\mathrm{Si}}}{a_{\mathrm{Si}}}\right],$$
(6)
where $`c_{11}`$ and $`c_{12}`$ are elastic constants for pure Si.

Strain produces shifts of the $`\mathrm{\Delta }`$ band proportional to the strain variables

$$\epsilon _{}(x)=\frac{a_{}(x)a_{\mathrm{Si}}}{a_{\mathrm{Si}}}\mathrm{and}\epsilon _{}(x)=\frac{a_{}(x)a_{\mathrm{Si}}}{a_{\mathrm{Si}}}.$$
(7)
with proportionality constants called the dilational and uniaxial deformation potentials, $`\mathrm{\Xi }_d^{(\mathrm{\Delta })}`$ and $`\mathrm{\Xi }_u^{(\mathrm{\Delta })}`$, respectively. Because of the anisotropic nature of the strain, the two $`\widehat{𝐳}`$ minima are shifted down relative to the $`\widehat{𝐱}`$ and $`\widehat{𝐲}`$ minima, resulting in a splitting of the $`\mathrm{\Delta }`$ conduction band. The net shifts with respect to the unstrained Si $`\mathrm{\Delta }`$ band are given by rieger93

$`\mathrm{\Delta }E^{(\mathrm{\Delta }_{})}(x)`$ $`=`$ $`(\mathrm{\Xi }_d^{(\mathrm{\Delta })}+{\displaystyle \frac{1}{3}}\mathrm{\Xi }_u^{(\mathrm{\Delta })})(2\epsilon _{}(x)+\epsilon _{}(x))+{\displaystyle \frac{2}{3}}\mathrm{\Xi }_u^{(\mathrm{\Delta })}(\epsilon _{}(x)\epsilon _{}(x)),`$ (8)
$`\mathrm{\Delta }E^{(\mathrm{\Delta }_{})}(x)`$ $`=`$ $`(\mathrm{\Xi }_d^{(\mathrm{\Delta })}+{\displaystyle \frac{1}{3}}\mathrm{\Xi }_u^{(\mathrm{\Delta })})(2\epsilon _{}(x)+\epsilon _{}(x)){\displaystyle \frac{1}{3}}\mathrm{\Xi }_u^{(\mathrm{\Delta })}(\epsilon _{}(x)\epsilon _{}(x)).`$ (9)
The first terms in Eqs. (8) and (9) are hydrostatic strain terms, which lower the conduction edge compared to unstrained Si. The second terms in Eqs. (8) and (9) produce the splitting, associated with uniaxial strain. To perform our calculations, we use the materials parameters given in Table 1. However we note that the deformation potentials, particularly $`\mathrm{\Xi }_d^{(\mathrm{\Delta })}`$, are very difficult to measure experimentally. Considerable disagreement exists in the literature as to the value and even the sign of $`\mathrm{\Xi }_d^{(\mathrm{\Delta })}`$ fischetti96 . The value given in Table 1 was reported (but not endorsed!) in Ref. fischetti96 , and provides energy band variations in general agreement with Refs. rieger93  and schaeffler97 . We arrive at the following strain-induced shifts of the conduction band edge for pure Si:

$$\mathrm{\Delta }E^{(\mathrm{\Delta }_{})}=0.86x\mathrm{and}\mathrm{\Delta }E^{(\mathrm{\Delta }_{})}=0.16x.$$
(10)
The corresponding shift in the relaxed barrier layers, due to the presence of Ge, was given in Eq. (4). Together, these results describe the conduction band offsets for the quantum well that are used in our simulations.

We now apply our results to two specific quantum well designs of interest for quantum computing. Design (1), shown in the inset of Fig. 1, is a version of that proposed by Vrijen et al. vrijen , in which electrons are trapped on donor ions (usually P), implanted in a semiconductor matrix. In that work, the quantum well is split into Ge-rich and Si-rich regions to facilitate single qubit operations. For simplicity, we consider here a uniform quantum well, formed of pure Si, with a single dopant ion located at the center of the well. In such a device, single qubit operations can be accomplished using a coded qubit scheme divincenzo00 . Design (2), proposed by Friesen et al. eriksson , is shown in Fig. 2. The confinement potential for the electrons is much softer than in design (1). Electrons are trapped vertically by the quantum well, and laterally by the electrostatic potential arising from lithographically patterned, metallic top-gates. Additionally, the quantum dot is tunnel-coupled to a degenerate doped back-gate. The dimensions for both designs are given in the figures.

The wave function of the bound electron is computed in the envelope function formalism, Eq. (3). Couplings between the different valleys are introduced through the perturbation theory described in Sec. IV. This procedure provides the specific values of $`\alpha _g^{(i)}`$ for the ground state, which we use in our calculations.

The abrupt conduction band offsets are handled by matching the ground state wave function $`\mathrm{\Phi }_g(\stackrel{}{r})`$ and $`_z\mathrm{\Phi }_g(\stackrel{}{r})`$ at the interfaces. (Remember that we have equated effective masses on both sides of interfaces.) Due to linear independence of the Bloch functions, the boundary conditions do not cause mixing of the envelope functions. Solutions of Eq. (3) and the analogous $`F_x`$ equation are obtained, using commercial 3D finite element software.

As will be seen in Sec. V, the key quantity for the computation of $`T_1`$ is $`f^2`$, which describes the probability for the bound electron to be on a Ge atom. Ge is associated with reduced coherence times, by virtue of its large spin-orbit coupling. Referring to Eq. (1), one deduces that this probability may be expressed as

$$f^2=x\left[4(\alpha _g^{(x)})^2_{\mathrm{\Omega }_b}d^3rF_x^2(𝐫)+2(\alpha _g^{(z)})^2_{\mathrm{\Omega }_b}d^3rF_z^2(𝐫)\right],$$
(11)
where $`\mathrm{\Omega }_b`$ is the volume outside the quantum well, if the well is pure Si. The subscript $`g`$ refers to the ground state. The term in square brackets reflects the probability of finding the bound electron in a barrier region, while $`x`$ gives the probability that the electron is on a Ge site.

Fig. 1 shows the results of our calculations for $`f^2`$ in design (1), as a function of the Ge concentration, $`x`$, in the Si1-xGex barriers. For $`x>0.02`$, $`f^2`$ decreases with $`x`$ for two reasons. First, as $`x`$ increases, the conduction band offset at the quantum well also increases, allowing less of the wavefunction to penetrate the barrier. Second, the spatial extent the electron in the $`\widehat{𝐳}`$ direction is greater for $`F_x`$ than $`F_z`$, because of the anisotropic effective mass . However, less of $`F_x`$ is mixed into the wavefunction for large $`x`$, since $`\alpha _x`$ becomes very small. For $`x0.02`$, $`f^2`$ drops quickly to zero, due to the absence of Ge in the barriers. In the actual design of Vrijen et al., vrijen  there is Ge in the active layer. To give an idea of the effect of this, we include an equivalent value of $`f^2`$ for such a structure.

Fig. 2 shows results of the $`f^2`$ calculation for design (2), as a function of the quantum well thickness, $`z`$. To perform the calculations, we have considered a fixed Ge concentration, $`x=0.05`$, and taken the limit of large strain, so that $`\alpha _g^{(\pm x)}\alpha _g^{(\pm y)}0`$ and $`\alpha _g^{(\pm z)}1/\sqrt{2}`$ for the ground state. As $`z`$ increases, less of the wavefunction penetrates the barrier regions, causing $`f^2`$ to decrease.

## III Spin relaxation due to coupling to phonons

In this section we give the method for calculating $`T_1`$, the spin-flip time of a spin qubit in the ground orbital state due to emission or absorption of a phonon, following the logic used by Hasegawa hasegawa  and Roth roth  for bulk Si.

Consider a single impurity with a unit positive charge, such as a phosphorus atom, at the origin. In the absence of central cell corrections, there is a 12-fold degenerate ground state, including spin. This valley degeneracy of the ground state is reduced to two by these corrections, and the splitting between the 2-fold spin-degenerate ground state and the higher states is of order $`\mathrm{\Delta }E10meV.`$ We shall discuss the detailed linear combinations ($`\alpha `$-values) of the states below, as the coefficients giving the various valley amplitudes play an important role in the calculation of matrix elements. These twelve states may all be thought of as hydrogenic $`1s`$ states. The splitting of $`1s`$ and $`2s`$ is about $`30meV`$, larger than the $`10meV`$ valley splittings. Let us now split the 2-fold degenerate ground state by applying a DC magnetic field in the $`z`$-direction. The transition rates between these states are denoted by $`W_{}`$ and $`W_{}`$. The relaxation time $`T_1`$ is defined by $`1/T_1=`$ $`W_{}+W_{}`$.

The transitions are caused by phonons, but there are important approximate symmetries that suppress these transitions. These are: (1) spin rotation symmetry, meaning that the electron spin cannot be flipped by a phonon; this symmetry is broken by SOC; (2) time-reversal symmetry, meaning that one state cannot be changed into its time-reversed partner by emission or absorption of a phonon; this symmetry is broken by the external magnetic field; (3) point group symmetries; these are partially broken when strain is applied.

The spin rotation symmetry would rule out phonon-mediated transitions between the two states entirely if there were no spin-orbit coupling (SOC). This means that the effects of SOC on the wavefunctions, even though these effects are small in relatively low-$`Z`$ Si, must be taken into account. When we refer to a state as $``$ or $``$, these symbols must be taken to refer to the majority spin content of the state, not to a pure spin state. Transition rates are roughly proportional to $`(g2)^2`$, \[more precisely $`(g_lg_t)^2`$, where $`g_t(g_l)`$ is the transverse (longitudinal) $`g`$-factor, see below for definitions\].

The time-reversal symmetry implies that transitions cannot take place directly between Kramers-degenerate states even in the presence of spin-orbit coupling. The direct phonon-mediated transitions between the two states of interest to us are strongly suppressed by this approximate symmetry. It is broken only by the external field $`H`$. The fastest processes then involve virtual excitation to higher-energy states that are mixed into the ground state by $`H`$. Hence $`1/T_1`$ involves a factor $`\left(\mu _BH/\mathrm{\Delta }E\right)^2`$. There is an additional factor of $`H^2`$ from the phonon density of states, giving an overall rate $`1/T_1H^4`$ in the limit of small $`H`$.

The point group symmetry is reduced from cubic to tetragonal under strain. This has complex effects that we will explain below.

Before giving actual calculations, we summarize those differences between the electrons in donor impurity states and in an artificial dot that affect $`T_1`$. The most obvious is the single-particle potential that binds the electron. The gate potential is much smoother than the hydrogenic potential of the impurity. This implies that the corrections to the effective mass approximation are much weaker, and $`\mathrm{\Delta }E`$ will be much reduced. It is difficult to compute the energy splittings precisely, but considerations based on the method of Sham and Nakayama sham79  give splittings in the range $`0.050.1meV`$ in the structure of Friesen et al. eriksson . This increases the relaxation rate. (In fact a naive estimate of the enhancement is a factor of 400.) On the other hand, the structures we consider have strong lattice strain. This partly lifts the valley degeneracy and also reduces the matrix elements, which decreases the rate. Another aspect of some of the proposed designs is the presence of Ge with its much stronger SOC. This will act to decrease the spin relaxation time.

## IV Pure Si quantum dots

We consider first the case of pure Si under uniaxial strain. The ingredients of the calculation are as follows.

From Sec. II we have the solutions to the Schrödinger equation $`[_0+V_g(\stackrel{}{r})]\mathrm{\Phi }_{n0}(\stackrel{}{r})=E_n\mathrm{\Phi }_{n0}(\stackrel{}{r})`$. $`_0`$ is the unperturbed crystal Hamiltonian without SOC and it has the full space group symmetry. $`V_g(\stackrel{}{r})`$ is the gate and/or impurity potential.

To calculate $`T_1`$, we must also include SOC, which we treat as a perturbation: $`H_{SOC}=\lambda _{Si}_\stackrel{}{R}\stackrel{}{L}_\stackrel{}{R}\stackrel{}{S}_\stackrel{}{R}`$. The resulting states $`\mathrm{\Phi }_n(\stackrel{}{r})`$ are two-fold degenerate because of time-reversal symmetry. Let us denote these states as $`\mathrm{\Phi }_n(\stackrel{}{r})`$ and $`\mathrm{\Phi }_n(\stackrel{}{r})`$. They are not eigenstates of spin, so the arrows denote pseudo-spin. We may define the pseudo-up state as that which evolves from the spin-up state as spin-orbit coupling is turned on adiabatically, and similarly for the pseudo-down state. Because of valley and pseudospin degeneracy, there are 2 ground states $`\mathrm{\Phi }_{gs}(\stackrel{}{r})`$ and 10 excited states $`\mathrm{\Phi }_{rs}(\stackrel{}{r}).`$ The 12-fold degeneracy in the effective mass approximation is broken by central-cell corrections in the impurity case and smaller corrections in the quantum dot.

There is also the Zeeman Hamiltonian of the external field $`_Z=\mu _B\stackrel{}{B}(\stackrel{}{L}+2\stackrel{}{S}).`$ In the field $`\mathrm{\Phi }_n(\stackrel{}{r})`$ and $`\mathrm{\Phi }_n(\stackrel{}{r})`$ are no longer degenerate. Note that the energy splitting may depend on the direction of $`\stackrel{}{B}`$.

Finally, we have the electron-phonon coupling Hamiltonian $`_{ep}.`$ A phonon represents a time-dependent perturbation. This will create transitions whose rate is given by the Fermi golden rule. We are interested in the transitions between $`\mathrm{\Phi }_g(\stackrel{}{r})`$ and $`\mathrm{\Phi }_g(\stackrel{}{r}).`$ However, $`\mathrm{\Phi }_g(\stackrel{}{r})|_{ep}|\mathrm{\Phi }_g(\stackrel{}{r})=0`$ in the absence of the external field. Thus we need to calculate in next order in perturbation theory using an effective Hamiltonian

$$^{}=\underset{rs}{}\frac{1}{E_gE_r}\left\{\left[_Z|\mathrm{\Phi }_{rs}(\stackrel{}{r})\mathrm{\Phi }_{rs}(\stackrel{}{r})|_{ep}\right]+\left[_{ep}|\mathrm{\Phi }_{rs}(\stackrel{}{r})\mathrm{\Phi }_{rs}(\stackrel{}{r})|_Z\right]\right\}$$
(12)
Here $`r`$ runs over the excited states, $`r=2,..,6`$. $`s=,`$.

The relaxation time is given by $`1/T_1=W_{}+W_{},`$ where $`W_{}=(2\pi /\mathrm{})\times _{\stackrel{}{q}\lambda }\left|\mathrm{\Phi }_g(\stackrel{}{r})|^{}|\mathrm{\Phi }_g(\stackrel{}{r})\right|^2\delta (E_gE_g\mathrm{}\omega _{\stackrel{}{q}\lambda })`$ $`[1+n(\omega _{\stackrel{}{q}\lambda })]`$ is the rate for transitions from the higher-energy (pseudo-spin down) state to the lower-energy (pseudo-spin up) state and $`W_{}=(2\pi /\mathrm{})_{\stackrel{}{q}\lambda }\left|\mathrm{\Phi }_g(\stackrel{}{r})|^{}|\mathrm{\Phi }_g(\stackrel{}{r})\right|^2\delta (E_gE_g\mathrm{}\omega _{\stackrel{}{q}\lambda })n(\omega _{\stackrel{}{q}\lambda })`$ is the rate for transitions from the lower-energy (pseudo-spin up) state to the higher-energy (pseudo-spin down) state where the sum is over phonon modes $`\stackrel{}{q}\lambda `$ with energies $`\omega _{\stackrel{}{q}\lambda }`$ . A thermodynamic average over the lattice states has been taken. It yields the Bose occupation factors $`n(\omega _{\stackrel{}{q}\lambda })`$ for the phonons.

The matrix elements of $`^{}`$ are computed as follows.

The expectation value of the external field $`H_Z`$ can be written as

$`\mathrm{\Phi }_{ns}(\stackrel{}{r})|H_Z|\mathrm{\Phi }_{n^{}s^{}}(\stackrel{}{r})`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{6}{}}}\alpha _n^{(i)}{\displaystyle \underset{j=1}{\overset{6}{}}}\alpha _n^{}^{(j)}\mu _B\stackrel{}{B}g^{(i)}\stackrel{}{\sigma }_{s,s^{}}\delta _{ij}`$ (13)
$`=`$ $`\mu _B\stackrel{}{B}\left[{\displaystyle \underset{i=1}{\overset{6}{}}}\alpha _n^{(i)}\alpha _n^{}^{(i)}g^{(i)}\right]\stackrel{}{\sigma }_{s,s^{}}`$
$``$ $`\mu _B\stackrel{}{B}𝐃_{nn^{}}\stackrel{}{\sigma }_{s,s^{}}`$
and the tensor $`g^{(i)}`$ is the effective $`g`$ factor at the $`i`$th valley. This equation defines the tensor $`𝐃_{nn^{}}`$ that characterizes the coupling of the various states by the external field. The principal axes of the $`g`$ tensor are the same as that of the effective mass tensor at the $`i`$th valley. It has the form

$$g^{(\pm x)}=\left(\begin{array}{ccc}g_l& 0& 0\\ 0& g_t& 0\\ 0& 0& g_t\end{array}\right),g^{(\pm y)}=\left(\begin{array}{ccc}g_t& 0& 0\\ 0& g_l& 0\\ 0& 0& g_t\end{array}\right),g^{(\pm z)}=\left(\begin{array}{ccc}g_t& 0& 0\\ 0& g_t& 0\\ 0& 0& g_l\end{array}\right).$$
(14)
There are only two independent constants.

A simple example of the diagonal part of the D tensor is that for the ground state of an impurity in the unstrained lattice when the central cell corrections are included. Then we have $`\alpha _g^{(i)}=1/\sqrt{6}`$ and $`𝐃_{gg}`$ is proportional to the unit matrix:

$$\mathrm{\Phi }_{gs}(\stackrel{}{r})|H_z|\mathrm{\Phi }_{gs^{}}(\stackrel{}{r})=g_g\mu _B\stackrel{}{B}\stackrel{}{\sigma }_{ss^{}},$$
(15)
with

$$g_g=\frac{2}{3}g_t+\frac{1}{3}g_l.$$
(16)
The matrix elements between ground and excited states have the form

$$\mathrm{\Phi }_{gs}(\stackrel{}{r})|H_Z|\mathrm{\Phi }_{rs^{}}(\stackrel{}{r})=g^{}\mu _B\stackrel{}{B}𝐃_{gr}\stackrel{}{\sigma }_{ss^{}},$$
(17)
with

$$g^{}=\frac{1}{3}(g_lg_t)$$
(18)
and the tensor $`𝐃_r`$ is defined by

$$𝐃_{gr}=3\underset{i=1}{\overset{6}{}}\alpha _g^{(i)}\alpha _r^{(i)}\widehat{k}^{(i)}\widehat{k}^{(i)}$$
(19)
where $`\widehat{k}^{(i)}`$ is the ”local” anisotropy axis. If the original $`g`$ were isotropic, then $`g^{}=0`$ and there would be no coupling between different states and no spin relaxation.

If the lattice is strained, then the $`\alpha `$ coefficients become strain-dependent and the general expression for $`𝐃`$ from Eq. 13 must be used. Uniaxial strain lifts the degeneracy of the valleys. We include this effect in the Hamiltonian and it determines the proper combinations of the $`\alpha _n^{(i)}`$ defined in Eq. 1. These then feed into $`𝐃_{gr}.`$ As a function of strain the $`\alpha _n^{(g)}`$ for the ground state cross over from the completely symmetric combination $`\alpha _g^{(i)}=1/\sqrt{6}`$ to the combination $`\alpha _g^{(\pm x)}=\alpha _{(g)}^{\pm y}=0`$, $`\alpha _{(g)}^{\pm z}=1/\sqrt{2}`$ in the limit of large strain.feher

The phonons involved are just the acoustic ones, one longitudinal and two transverse - these are the only ones with low enough energy to play a role in relaxing the spins. The matrix elements of the electron-phonon interaction are only nonzero within one valley and for a single phonon mode they are conventionally parametrized as:

$$\mathrm{\Psi }_{c\stackrel{}{k}s}^{(i)}(\stackrel{}{r})|H_{ep}^{(\stackrel{}{q}\lambda )}|\mathrm{\Psi }_{c\stackrel{}{k}^{}s}^{(i)}(\stackrel{}{r})=ib_{\stackrel{}{q}\lambda }\widehat{e}_\lambda (\stackrel{}{q})(\mathrm{\Xi }_d1+\mathrm{\Xi }_u\widehat{k}^{(i)}\widehat{k}^{(i)})\stackrel{}{q}+h.c.$$
(20)
near the $`i`$th valley, where $`\stackrel{}{q}=\stackrel{}{k}\stackrel{}{k}^{}`$ and $`\widehat{e}_\lambda `$ is the polarization vector. $`b_{\stackrel{}{q}\lambda }`$ destroys a phonon with wavevector $`\stackrel{}{q}`$ and polarization $`\lambda `$. Once again, we see that the interaction can be characterized by just 2 parameters, in this case $`\mathrm{\Xi }_d`$ and $`\mathrm{\Xi }_u`$, as already defined in Eq. 7. Performing the integration over the envelope function at wavector $`\stackrel{}{q}`$ now gives

$`\mathrm{\Phi }_{gs}(\stackrel{}{r})|H_{ep}^{(\stackrel{}{q}\lambda )}|\mathrm{\Phi }_{gs^{}}(\stackrel{}{r})_\stackrel{}{q}`$ $`=`$ $`(\mathrm{\Xi }_d+{\displaystyle \frac{1}{3}}\mathrm{\Xi }_u)A(\stackrel{}{q})\left(b_{\stackrel{}{q}\lambda }+b_{\stackrel{}{q}\lambda }^{}\right)\delta _{s,s^{}}`$ (21)
$`\mathrm{\Phi }_{gs}(\stackrel{}{r})|H_{ep}^{(\stackrel{}{q}\lambda )}|\mathrm{\Phi }_{rs^{}}(\stackrel{}{r})_\stackrel{}{q}`$ $`=`$ $`{\displaystyle \frac{1}{3}}\mathrm{\Xi }_uA(\stackrel{}{q})\left(i\widehat{e}_\lambda (\stackrel{}{q})𝐃_{gr}\stackrel{}{q}b_{\stackrel{}{q}\lambda }i\widehat{e}_\lambda ^{}(\stackrel{}{q})𝐃_{gr}\stackrel{}{q}b_{\stackrel{}{q}\lambda }^{}\right)\delta _{s,s^{}}`$ (22)
where

$$A^{(i)}(\stackrel{}{q})=\underset{\stackrel{}{k}}{}F^{(i)}(\stackrel{}{k}+\stackrel{}{q})F^{(i)}(\stackrel{}{k})=d^3rF^2(r)e^{i\stackrel{}{q}\stackrel{}{r}}.$$
(23)
Thus the electron-phonon interaction involves a form factor for the bound states. Since $`F`$ is normalized, we have $`A^{(i)}(\stackrel{}{q})1`$ when the wavelength of the phonon is much longer than the spatial extent $`a^{}`$ of the bound state: $`qa^{}<<1`$. The calculations of II indicate that this is the case. It is also independent of $`(i).`$

In the golden rule calculation, the energy denoninator $`\left(E_gE_r\right)^2`$ will suppress contributions from the excited states of $`V_g`$. Thus we will keep only states that are split off from the ground state by corrections to the effective mass approximation. This approach works very well in Si:P and should be even better for the quantum dot.

This produces the golden-rule transition rate

$`W_{}`$ $`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{}}}\left[{\displaystyle \frac{1}{3}}\mathrm{\Xi }_ug^{}\mu _B\right]^2{\displaystyle \underset{\stackrel{}{q}\lambda }{}}A^2(\stackrel{}{q})\delta (E_gE_g\mathrm{}\omega _{\stackrel{}{q}\lambda })a_{\stackrel{}{q}\lambda }a_{\stackrel{}{q}\lambda }^{}\times `$ (24)
$`\left|{\displaystyle \underset{r=2}{\overset{6}{}}}{\displaystyle \frac{\stackrel{}{B}𝐃_{gr}\stackrel{}{\sigma }_{}\widehat{e}_\lambda (\stackrel{}{q})𝐃_{gr}\stackrel{}{q}\delta _{ss^{}}}{E_gE_r}}\right|^2.`$
We approximate the phonon dispersion as $`\omega _{\stackrel{}{q}\lambda }=v_\lambda q`$. Setting $`A=1`$, performing the integral over the magnitude of $`\stackrel{}{q}`$, and repeating the calculation for $`W_{}`$, we obtain a total spin relaxation time

$`{\displaystyle \frac{1}{T_s}}`$ $`=`$ $`W_{}+W_{}={\displaystyle \frac{1}{8\pi ^2\rho \mathrm{}^4}}\left({\displaystyle \frac{g^{}\mu _BB\mathrm{\Xi }_u}{3}}\right)^2[2n(g\mu _BB)+1]g_g^3\mu _B^3B^3\times `$ (25)
$`{\displaystyle \underset{\lambda =1}{\overset{3}{}}}{\displaystyle \frac{1}{v_\lambda ^5}}{\displaystyle _0^{2\pi }}𝑑\varphi ^{}{\displaystyle _0^\pi }\mathrm{sin}\theta ^{}d\theta ^{}\left|{\displaystyle \underset{r=2}{\overset{6}{}}}\left[\stackrel{}{B}(\theta ,\varphi )𝐃_{gr}\stackrel{}{\sigma }_{}\right]\left[{\displaystyle \frac{\widehat{e}_\lambda (\stackrel{}{q})𝐃_{gr}\widehat{q}}{E_gE_r}}(\theta ^{},\varphi ^{})\right]\right|^2.`$
Here $`(\theta ,\varphi )`$ are the polar axes of the direction of $`\stackrel{}{B}`$ measured from the $`[100]`$ direction of the crystal. $`\rho `$ is the mass density. This is our basic result, in a form very similar to that given by Hasegawa hasegawa .

In Fig. 3 the effect of strain on the lowest energy levels is shown. Here the zero-strain splittings are those of a P impurity. A similar plot was given in Ref. feher , but at that time the correct zero-strain splittings were not known. In the dot, the energy splittings for the unstrained case would be two to three orders of magnitude smaller (see below for a discussion).

In Fig. 4, we show the effect of strain on the relaxation time of an electron in a P impurity potential. There are two effects: the overall increase of the energy denominators and the change of the ground state to a less symmetric valley weighting. This leads to a non-monotonic dependence of $`T_1`$ on strain, but in the region of interest, the effect of strain is to greatly increase $`T_1`$, since most proposed structures have $`s<1.`$ At large strain, only one energy denominator remains small, that between the ground state, symmetric in the $`\pm z`$ valleys, and the first excited state, antisymmetric in the $`\pm z`$ valleys. The overlap matrix D is very small between these two state. This reduction of the matrix element is the dominant effect.

Also of interest is the dependence of $`T_1`$ on the angle of the external field, since this may serve as a diagnostic tool in experiments to verify that the relaxation process is really due to spin-phonon coupling. As a function of strain, this dependence becomes highly anisotropic, as seen in Fig. 5. This is due to the elimination of all but the $`\pm z`$ valleys from the problem at high strain. We note that the limiting value $`T_1\mathrm{}`$ when the external field is along a crystal axis is cut off by intervalley scattering effects not included in the present calculation feher .

The change in the confining potential reduces the corrections to the effective mass approximation, as discussed in Sec. III which in turn reduces the energy denominators, leading to a decrease in $`T_1`$. In the Friesen et al. structure, we have estimated the splitttings $`E_rE_g`$ using the method of Sham and Nakayama sham79 , and they range from 0.05 to 0.1 meV, depending on the gate voltages that produce the potential. This may be incorporated into the calculation of the $`\alpha _n^{(i)}`$ by introducing a variable coupling $`\mathrm{\Delta }_c`$ that mixes the $`\pm x`$ and $`\pm y`$ valley wavefunctions with the $`\pm z`$ valley wavefunctions. For a precise definition of $`\mathrm{\Delta }_c`$, see Ref. feher . Approximately, however, $`\mathrm{\Delta }_c0.2(E_rE_g)<10^4eV`$. The results of the calculations are shown in Fig. 6. The graphs show the absolutely crucial role that strain plays in the determination of $`T_1`$. Because of the small energy denominators in the dot, the rate is extremely fast at small $`\mathrm{\Delta }_c`$ for the unstrained case: $`1/T_1\mathrm{\Delta }_c^2.`$

We also point out that in stuctures where $`T_1`$ appears to be too long for efficient preparation of the spins (as discussed in Sec. I) the deficiency can be made up by increasing the temperature $`T`$. As seen in Fig. 7, $`T_1`$ decreases very rapidly as $`T`$ is increased. Actually, this calculation even underestimates the decrease, since multiphonon processes begin to contribute at about 3K. However, it must be borne in mind that the temperature must be small enough that the initial state is essentially pure (all spins up, for example). A more practical method of decreasing $`T_1`$ if desired, would be to increase the magnetic field.

## V Structures containing Ge

The effect of alloying with Ge is to increase the SOC and hence to increase $`\left|g2\right|`$. First principles calculations of Ge impurities in Si have shown that there is little effect on the states near the bottom of the conduction band, though this is not necessarily the case for higher energy states in the band jaros . This is in accord with the isoelectronic character of the atoms in the alloy. It is then reasonable to employ the virtual crystal approximation (VCA). The approximation should be quantitatively accurate for small, (say $`<10\%`$) concentrations of Ge, but may be taken as a good qualitative guide also to higher concentrations. The Bloch functions in the absence of SOC satisfy

$$H_0\mathrm{\Psi }_{n\stackrel{}{k}s}=E_{n\stackrel{}{k}}\mathrm{\Psi }_{n\stackrel{}{k}s},$$
(26)
where

$$\mathrm{\Psi }_{n\stackrel{}{k}s}=\frac{1}{\sqrt{N_c}}\underset{j=1}{\overset{N_c}{}}\underset{l=1}{\overset{2}{}}\mathrm{exp}(i\stackrel{}{k}r_{jl})\underset{m=0}{\overset{3}{}}a_{lm}(n\stackrel{}{k})\varphi _{ms}(\stackrel{}{r}r_{jl})$$
(27)
where $`N_c`$ is the number of unit cells, $`j`$ labels the unit cells, $`l`$ labels the two positions in the unit cell, $`m`$ labels the 4 atomic states, $`\varphi _{ms}`$ are the atomic orbitals, and the coefficients $`a_{lm}(n\stackrel{}{k})`$ give the proper linear combination of atomic orbitals (LCAO’s) for the state at momentum $`\stackrel{}{k}`$ and band $`n`$. The normalization condition is

$$\underset{lm}{}\left|a_{lm}(n\stackrel{}{k})\right|^2=1.$$
(28)
In the VCA, the $`a_{lm}(n\stackrel{}{k})`$ are the same for the Si and Ge sites. The SOC Hamiltonian for the alloy is written as

$$H_{SOC}=\underset{\stackrel{}{R}}{}\lambda _\stackrel{}{R}\stackrel{}{L}_\stackrel{}{R}\stackrel{}{S}_\stackrel{}{R},$$
(29)
where $`\lambda _\stackrel{}{R}`$ is the SOC strength for Si (Ge) when $`\stackrel{}{R}`$ is a Si (Ge) site. The energy shift of an electron in the external field $`\stackrel{}{B}=(B_x,B_y,B_z)`$ can be written as

$$\mathrm{\Delta }E_{n\stackrel{}{k}s}(\stackrel{}{B})=\mu _B\underset{i=x,y,z}{}g_i(n\stackrel{}{k})B_is,$$
(30)
where

$$g_i(n\stackrel{}{k})=22\mathrm{}^2\left(\underset{\stackrel{}{R}}{}\lambda _\stackrel{}{R}/N_c\right)\underset{n^{}n,l,m,m^{}}{}\left|a_{lm}(n\stackrel{}{k})\right|^2\left|a_{lm^{}}(n^{}\stackrel{}{k})\right|^2\left(E_{n\stackrel{}{k}}E_{n^{}\stackrel{}{k}}\right)^1(\epsilon _{imm^{}})^2$$

and $`\epsilon _{imm^{}}`$ is the completely antisymmetric symbol with $`\epsilon _{123}=1`$.

Our approximation now consists in regarding the dependence on the Ge concentration as coming only in the term $`\left(_\stackrel{}{R}\lambda _\stackrel{}{R}/N_c\right)`$. Hence any component of the $`g`$-tensor is proportional to a weighted average of $`\lambda _{Si}`$ and $`\lambda _{Ge}.`$ In particular, consider the conduction band minimum at $`\stackrel{}{k}=(0,0,k_0)`$ in the compound Si1-xGe$`_x.`$ Then we have $`g_z(x)2=[g_z(0)2]\left[(1x)+x(\lambda _{Ge}/\lambda _{Si})\right],`$ where $`g_z(0)`$ is the value for pure silicon. The numbers, known from atomic physics, are $`\lambda _{Ge}/\lambda _{Si}=0.295/0.044=6.71`$, so $`g_z(x)2=\left[g_z(0)2\right]\left[(1x)+6.71x\right].`$ Similarly, we have $`g_y(x)2=\left[g_y(0)2\right]\left[(1x)+6.71x\right]`$ , etc. Only the overall scale of $`g`$ changes, not the anisotropy in the VCA.

The presence of Ge in the lattice breaks the translation symmetry, an effect that is neglected in the VCA. Taking this into account would lead to a $`g`$-factor that is averaged over momentum space rather than one that is evaluated only at the Si conduction band minima. If we take the L-minimum of pure Ge as an example, we would expect that the corrections to $`\left|g2\right|`$ would be greater if there is some averaging over momentum space. Thus a calculation going beyond the VCA would give corrections that would reduce $`T_1`$.

We show the results for differing Ge concentrations in the active region of the well in Fig. 8. In the Ge-poor regime $`x<0.85`$, the relaxation rate $`1/T_1`$ decreases fairly slowly. Thus designs with a Si-rich well will not suffer from very fast spin relaxation on the minority Ge sites.

By contrast, if we cross the critical point around $`x=0.85`$ where the conduction band minima switch positions, then there is a very rapid decrease of $`T_1`$. This is due to the much greater g-factor anisotropy of the  minima: $`g^{}0.4`$, (as compared to $`g^{}<0.001`$ in Si). The results for the Ge-rich region are much simpler than in the Si-rich region. Uniaxial strain does not affect the relative valley energies or the matrix elements parametrized by D. The only effect of going to the dot from the impurity is to reduce corrections to the effective mass approximation, which will strongly decrease the relaxation time. Indeed, if we take the impurity result for Ge:P from Ref. hasegawa , which is $`T_1=2.3\times 10^3s`$ and estimate the decrease in $`\mathrm{\Delta }E`$ as about a factor of 100, then we obtain $`T_110^6`$ to $`10^7s,`$ which is vastly shorter than that of Si-rich structures. The physics is the same as that which governs the divergence of the rate in Fig. 6 for the unstrained case,namely the much smaller energy denominators.

To make contact with Sec. II note that if there is Ge in the barrier regions, as is the case in most SiGe designs, then we will have an envelope function that weights Si and Ge differently, then $`x`$ in the VCA formulas must be replaced by $`f^2=_\stackrel{}{R}|f_\stackrel{}{R}|^2`$, where $`|f_\stackrel{}{R}|^2`$ is the amplitude that the electron is on a Ge site and the sum runs only over Ge sites, as described in Sec. II. This is the product of the Ge concentration in the layer, times the probability that the electron is on the layer, summed over layers. Referring to Figs. 1 and 2, we see that $`f^2`$ is always less than $`10^3`$ for structures that have pure Si active layers and Ge in the barrier layers. The effect on $`T_1`$ is small so these structures are ”phonon-safe”.

## VI Conclusion

Decoherence due to spin relaxation by emission and absorption of single phonons does not pose a major obstacle to quantum computation in quantum dot SiGe heterostructures using spin qubits, if these qubits are properly designed. This spin relaxation mechanism is the dominant one in the case of electrons in donor bound states, such as in Si:P.

In order that this source of decoherence be kept under control, however, certain conditions must be satisfied in the design of the structure. Unstrained SiGe alloys generally have relatively short spin relaxation times due to valley degeneracy. This degeneracy produces strong spin mixing in the Kramers-degenerate ground states of the quantum dot once the field is turned on. This in turn increases the electron-phonon matrix element that connects the two states and decreases the relaxation time. This effect, fundamentally due to spin-orbit coupling, is stronger in the quantum dot than in the impurity state. The corrections to effective-mass theory, stronger for the impurity potential, tend to lift the degeneracy.

The crucial role that strain plays, in Si-rich structures, is to lift the valley degeneracy. In this respect compressive uniaxial strain is the best, since the residual degeneracy is reduced to two valleys. The wavefunction for the ground state is then symmetric in the valley index, and the first excited state, the only one with a small energy denominator, is antisymmetric in this index. This fact strongly suppresses the spin mixing effect of the external field. Thus a workable design should include uniaxial strain that exceeds a critical value. This value is small enough that there will be sufficient strain in most practical structures in which the active layer is Si rich and sandwiched between relatively Ge-rich layers. The active layer must be thin enough to avoid dislocation-mediated relaxation.

Spin-orbit coupling is stronger in Ge than in Si due to the higher atomic number. This suggests that Si-rich structures are to be preferred, and the calculations bear this out. The relaxation time decreases with increasing $`x`$ in a Si1-xGex alloy. In fact, if $`x`$ increases beyond $`x=0.85`$, decoherence becomes very strong, and quantum dots in these Ge-rich structures, such as those proposed by Vrijen et al. vrijen , may well run into difficulties for this reason. A key point for these structures is that uniaxial stress does not lift the valley degeneracy, since the valleys have moved off the crystal axes. Strain-induced suppression of spin relaxation cannot occur.

Beyond design issues, fabrication quality is also important. Even a small concentration of magnetic impurities will negatively impact the spin relaxation. In modern semiconductor technology, however, concentrations of magnetic impurities much less than 0.1 ppm are routine. The impact of lattice imperfections is less clear. These will generally act to lower the symmetry of the system, and to lessen the accuracy of the effective-mass approximation. As we have seen in the context of the impurity calculations, these effects generally increase the spin relaxation time because of reduced degeneracy. On the other hand, localized bound states can form at such imperfections. If this results in free electron spins, it could have seriously negative effects on the relaxation time. Similar effects would result from lattice defects that produce localized phonon modes.

Finally, it appears that quantum dot designs will need rather low temperatures in order to operate. The impurity experiments clearly indicate that multi-phonon relaxation increase rapidly with temperature tgcastner , with a crossover to this regime at about 3 K. We speculate that this crossover temperature is pushed up when the effects of strain are included, but we have not yet done the requisite calculations.

We would like to acknowledge useful discussions with our collaborators at the University of Wisconsin-Madison and with Xuedong Hu. This work was supported by the National Science Foundation under the Materials Theory Program, Grant No. DMR-0081039, the Materials Research Science and Engineering Center Program, Grant No. DMR-0079983, and the Army Research Office.


text:
# Microstructure and pinning properties of hexagonal-disc shaped single crystalline MgB2

New address: \]Tokura Spin Superstructure Project, ERATO, JST, AIST, Tsukuba Central 4 Tsukuba, Ibaraki 305-8562, Japan Url: \]http://unit.aist.go.jp/cerc/ Url: \]http://www-psc.postech.ac.kr

## Abstract

We synthesized hexagonal-disc-shaped MgB2 single crystals under high-pressure conditions and analyzed the microstructure and pinning properties. The lattice constants and the Laue pattern of the crystals from X-ray micro-diffraction showed the crystal symmetry of MgB2. A thorough crystallographic mapping within a single crystal showed that the edge and c-axis of hexagonal-disc shape exactly matched the (10-10) and the (0001) directions of the MgB2 phase. Thus, these well-shaped single crystals may be the best candidates for studying the direction dependences of the physical properties. The magnetization curve and the magnetic hysteresis for these single crystals showed the existence of a wide reversible region and weak pinning properties, which supported our single crystals being very clean.

MgB2  preprint: PRB-MgB2 Single crystal

The recent discoveryfind  of superconductivity in MgB2 has attracted great scientificKangHall ; BudkoPRB  and industrial Kangfilm ; Larbal  interest. Even though basic issues such as the carrier type KangHall  were addressed immediately, conflicting reports still exist, especially on the transport properties. For example, the residual resistivity ratio (RRR) of bulk MgB2 ranges from 1.2 to 30. JungNKim

Higher quality bulk MgB2 (called highRRR-MgB2) has been claimed to have a higher values of RRR $`(2025)`$, a low residual resistivity ($`\rho (40`$ K$`)<1`$ $`\mu \mathrm{\Omega }`$cm), a higher magneto-resistance (MR), and a resistivity upturn at low temperature under high magnetic field. BudkoPRB  Insulating impurities and/or local strains have been suggested as possible origins for these different observations.Zhu ; Chu1 ; Chu2 ; Chu3 ; 0598  However, very recently, the existence of unreacted Mg successfully explained the unusual enhancement of the RRR and the MR in polycrystalline MgB2.JungNKim

Compared to the wide distribution of the RRR for polycrystals, the RRR in the ab-plane resistivity for single crystals has a narrow distribution.Lee ; Xu  The RRR of single crystals ranges from 5 to 6, which is much smaller than the values for highRRR-MgB2 bulk samples, but similar to the values for well-prepared polycrystals.JungNKim ; Nesterenko1 ; Nesterenko2 ; Fuchs  The anisotropy factor $`(\gamma 3)`$ and the Debye temperature $`(\mathrm{\Theta }1160\pm 60`$ K) are also consistent among various reports.MR

It is interesting to note that the magnetic properties are somewhat different for different single crystals.Lee ; Xu ; Lima  The superconducting transition width in the zero-field-cooled magnetization for single crystals as large as a few hundred $`\mu `$m is a little bit broader than expected. The ratio of the low-field magnetization in the field-cooled (FC) state to that in the zero-field-cooled (ZFC) state gives a rough indication of the pinning strength, and this ratio is quite different for different single crystals. The magnetic hysteresis for single crystals having very weak pinning can be used to probe impurities. In one report, the magnetic hysteresis contained a significant amount of paramagnetic component. In another report, the magnetic hysteresis for an aggregation of single crystals showed a large ferromagnetic contribution and a significant irreversible magnetization.Lima

For the above reasons, it is necessary to confirm the quality of single crystals in more detail before performing the main measurements. The existence of impurities and structural imperfections on a microscopic scale can result in diverse physical properties. Here, we report the growth of, as well as X-ray micro-diffraction and magnetization measurements for, MgB2 single crystals with hexagonal-disc shapes and shiny surfaces. Our single crystals are unique as far as the shape is concerned. The diagonal length and the thickness for the largest crystal was about 100 $`\mu `$m and 10 $`\mu `$m, respectively. The crystallinity was thoroughly identified by using the Laue pattern in the X-ray micro-diffraction measurement. Both the edge and the c-axis of the hexagonally shape disc were found to match the crystal symmetry. The magnetization study showed that pinning was very weak for our hexagonal-disc-shaped single crystals.

Two different procedures were used to grow the single crystals, and in both cases, excess Mg was critical for the growth of single crystals. The first involved a two-step method in which already synthesized pieces of MgB2 bulkJungHPsynthesis  were used as a seed material. They were heat treated in a Mg flux inside a Nb tube, which was sealed in an inert gas atmosphere. Then, the Nb tube was put inside a quartz tube, which was sealed in vacuum. The quartz tube was heated for one hour at 1050C, cooled very slowly to 700C for five to fifteen days, and then quenched to room temperature. The temperature dependences of the ab-plane resistivities for $`H//ab`$ and $`H//c`$ and high-resolution transmission electron microscope images have already been reported for the single crystals synthesized using this procedure.SingleKim

In the one-step method, 1 : 1 mixtures of Mg and B powders were well ground and pressed into a pellet. Then, the pellet was placed in a Ta capsule. This capsule was put in a high-pressure cell equipped with a graphite heater. Heat treatment was done inside a 14-mm cubic multi-anvil-type press under 3 GPa.JungHPsynthesis  The heating temperature was around $`1300^{}`$C. The temperature was maintained for about 30 minutes, and then was slowly lowered to $`800900^{}`$C. The final product was a pellet containing a mixture of single-crystalline MgB2 and Mg flux.

It was found that the one-step method usually gave larger crystals, which was quite advantageous for the magnetization study in this research. The crystal images were observed using a polarizing optical microscope and a field-emission scanning electron microscope (SEM). We successfully separated single-crystalline MgB2 from the Mg flux by using a thermo-mechanical spinning method. This method is possible due to the fact that the melting (and/or decomposition) temperature of MgB2 is higher than that of Mg.

Single crystals with sizes of tens of $`\mu `$m were selectively handled by using a homemade micro-tweezers and were fixed on Si substrates by using a photoresist as an epoxy. Since the volume of one crystal made in the one-step process was still rather small, we gathered about 200 single crystals on one substrate with their c-axes aligned perpendicular to the substrate surface.Lima  For the X-ray micro-diffraction measurements, several crystals were fixed at the center of Cu crosshairs on the substrate, as shown in Fig. 1(d).SingleKim  The Cu crosshairs was used for the only purpose to facilitate the location of the sample by looking at the Cu fluorescence. The magnetic properties were measured with a SQUID magnetometer (Quantum Design, MPMS-XL).

The instrument used at the Advanced Light Source (ALS) for X-ray micro-diffraction is capable of producing a submicron-size X-ray microbeam and with submicron spatial resolution can probe the local texture in a single crystal.NIMA  The sample was positioned using the Cu fluorescence signal detected from the Cu crosshairs on the Si substrate by using a high-purity Ge ORTEC solid-state detector connected to a multichannel analyzer. The crystal orientation with respect to the substrate can be determined with an accuracy of 0.01 degree.

Figure 1(a) shows a typical SEM image for a MgB2 single crystal synthesized in a two-step process. The crystals observed by using a polarizing optical microscope had hexagonal-disc shapes with edge angles of 120 degrees and very flat and shiny surfaces. The sizes of the crystals were about $`2060`$ $`\mu `$m in diagonal length and $`26`$ $`\mu `$m in thickness. Figure 1(b) shows a magnified view of the upper corner of another crystal. The smooth surfaces and the sharp edges confirm that our small crystals had a very low probability of having mosaic aggregates of nanocrystals either along the $`ab`$-plane or along the $`c`$-axis; thus, we had a better chance to study their intrinsic properties. The single crystals reported so far, except for those reported by us, have had irregular shapes. Lee ; Xu ; Lima

Figure 1(c) shows a polarizing optical microscope image of a MgB2 single crystal synthesized in a one-step process, which resulted in larger crystals. The sizes of these crystals were $`30100`$ $`\mu `$m in diagonal length, so they were picked for measurements of the magnetic properties and the X-ray micro-diffraction.

The crystal structure was identified by using white beam X-ray micro-diffraction measurements. After positioning these single crystals on the substrate in Fig. 1(d), a $`100`$ $`\mu `$m $`\times 100`$ $`\mu `$m region between the Cu crosshairs was scanned with a step size of 2 $`\mu `$m. At each step, the Laue pattern (together with the Cu K fluorescence signal ) was collected with a BRUKER 6000 CCD camera which has an active area of $`9\times 9`$ cm and was placed about 4 cm above the sample. \[2500 images, 1024 pixel $`\times 1024`$ pixel mode\] The exposure time at each step, was 1 second. An example of a Laue pattern obtained from a MgB2 single crystal is shown in Fig. 2. The Laue patterns are consistent with a hexagonal MgB2 structure ($`a=3.086`$ Å, $`c=3.524`$ Å, Space Group number =191, Ref. 1). The silicon reflections from the substrate were digitally removed. Typically, more than 20 reflections with energies ranging from 5 to 14 keV were indexed successfully. The (0005) reflection in the center of the pattern in Fig. 2 corresponds to the direction of the normal to the crystal surface. This confirms that the surface plane normal is along the c-axis. Moreover, the hexagonal edges of the crystals were found to match the $`<`$1,0,-1,0$`>`$directions within a fraction of a degree resolution. Thus, the shapes of the crystals in the microscope image followed the MgB2 crystal symmetry, which will be quite useful for any research of the direction dependencies of the physical properties in MgB2.

Indexing the Laue patterns in Fig. 2(a) allowed us to calculate the complete orientation matrix of the X-ray illuminated volume. A finer step size of 1 $`\mu `$m was used for the white-beam scan. The orientation variations inside the single crystal shown in the right bottom corner of Fig 1(d) are shown in Figs. 2(b) and 2(c). Figure 2(b) is the out-of-plane orientation variation calculated as the angle between the c-axis and the normal to the surface of the silicon substrate. The out-of-plane variation was about 0.2 degrees between the light orange and the red regions. The out-of-plane orientation shows a variation of about 0.2 degrees between the bottom and the top parts, indicating a slight bending of the crystal (which might be due to the photoresist used as an epoxy). Figure 2(c) shows the in-plane orientation variation calculated as the angle between the measured a-axis (or b-axis) and a reference directions. The variation was about 0.4 degrees between the light blue-green and the green regions. The in-plane orientation also showed some inhomogeneities of up to about 0.2 degrees.

These results demonstrate that the orientation of crystal axis of our hexagonal-disc-shaped single crystals was perfect, within 0.2 degrees. A recent study showed that (0001) twist grain-boundaries, formed by rotations along the $`c`$-axis (typically by about 4 degrees), were the major grain boundaries in polycrystalline MgB2.Zhu  This kind of grain boundary was attributed to the weaker Mg-B bonding.Zhu

To study the bulk nature of the superconductivity, we measured the magnetization curve $`M`$($`T`$) and the magnetic hysteresis $`M`$($`H`$). Figure 3(a) shows the magnetization curve $`M`$($`T`$) measured at 20 Oe in the ZFC and the FC states. The $`T_c`$ onset was about $`38`$ K. The FC signals for fields parallel and perpendicular to the c-axis were larger than 60% of the ZFC signals, which suggested that the pinning was very weak. The overall shape was not much different for these two orientations. The different values of the magnetic moment for the different field directions were due to a demagnetization effect caused by the planar-disc shapes of the crystals.SingleKim  The transition width in the $`M(T)`$ curve was much narrower than that for single crystals made in the two-step method.SingleKim  The volumes of the crystals made in the one-step process were about one order of magnitude larger than the volumes of the crystals made in the two-step process.

The weak pinning was demonstrated by the magnetization at 5000 Oe and the magnetic hysteresis $`M`$($`H`$) at 5 K, as shown in Fig. 3(b) and (c), respectively. With $`(M(H^{})M(H^+))=0.1`$ emu/cm3 as the criterion for the reversible point, the values of $`T_{\mathrm{irr}}(H=5000`$ Oe) were 26 and 24 K for $`H//c`$ and $`H//ab`$, respectively. At each temperature, the upper critical fields obtained from the resistivity measurements were about 1.7 and 6 Tesla for our well-shaped single crystals.SingleKim  Thus, a very wide reversible region existed in our single crystals, especially for $`H//ab`$. \[$`H_{\mathrm{irr}}`$ for bulk MgB2 at $`T=25`$ K was about 3 TFuchs , which is about 6 times higher than that of single crystalline MgB2. \]

The weak pinning shown by the wide reversible region was also consistent with the results of the $`M(H)`$ measurement. The $`M(H)`$ at 5 K in Fig. 3(c) showed negligible paramagnetic or ferromagnetic background up to 5 Tesla. Reversible magnetization was dominant for the $`M(H)`$ of our single crystals, which was quite consistent with fact that the pinning in our single crystals was weak and with the existence of wide reversible region demonstrated by the lower $`H_{\mathrm{irr}}`$ value in Fig 3(b). The different slopes of the $`M(H)`$ curves at the starting low fields ($`H<300`$ Oe) for the different field directions were due to different demagnetization factors.

Previously, the $`M(H)`$ for irregular shaped single crystals was shown to have a non-negligible irreversible contributionLima  and a quite large paramagnetic and/or ferromagnetic background signal, which might have been due to the existence of impurities.Xu ; Lima  The results in Fig. 3 indicate that the strong bulk pinning previously reported for polycrystalline kim  and thin filmsKangfilm  might be due to entirely extrinsic pinning sites, such as grain boundaries and crystallographic defects.Pissas  This is consistent with the absence of core pinning, even at $`T0.5\times T_c`$, for bulk sample.Larbal

In summary, we report the structural and magnetic properties of MgB2 single crystals with hexagonal-disc shapes. The X-ray micro-diffraction showed that the hexagonal-disc shape of the single crystal followed the crystallographic symmetry, which will be very useful for studying orientation-dependent physical properties. The magnetization curve and the magnetic hysteresis provided consistent evidence that our single crystals were very clean and had very weak pinning: the large FC magnetization and the sharp transition of the ZFC magnetization at low field $`M(T)`$, the wide reversible region, and the $`M`$($`H`$, $`T=5`$ K) dominated by the reversible magnetization without any significant paramagnetic and/or ferromagnetic contribution.

###### Acknowledgements.

We appreciate M. H. Kim and S. R. Jung for their help in handling of tiny crystals. This work is supported by the Ministry of Science and Technology of Korea through the Creative Research Initiative Program. The Advanced Light Source is supported by the Director, Office of Science, Office of Basic Energy Sciences, Materials Sciences Division, of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098 at Lawrence Berkeley National Laboratory.


text:
# Molecular simulation and theory of a liquid crystalline disclination core

## I Introduction

Liquid crystals are at the heart of a range of technological devices, which rely on the ability to manipulate the direction of preferred molecular alignment, the *director* $`𝐧`$, with electric and magnetic fields, and through coupling to surfaces. The practical performance of such devices, and the theoretical description of director distortions, both rely on smooth variation of the director with position in space; this is well accounted for by the Frank free energy frank.fc:1958.a ; degennes.pg:1995.a  which is an expansion in squared gradients of the director field $`𝐧(𝐫)`$, parametrized by the splay ($`K_{11}`$), twist ($`K_{22}`$) and bend ($`K_{33}`$) elastic constants. The director $`𝐧`$ is the principal axis of the local second-rank ordering tensor $`𝐐`$, which characterizes the nematic state; such a description is intrinsically uniaxial, and neglects variation of the *degree* of ordering with position, of the kind which occurs near bounding surfaces and around topological defects. Such defects are treated as singularities in the director field; for a better description, it is necessary to replace the director field by one which allows some variation of the relevant order tensor components over relatively short length scales. The relevant region near the defect is called the *core*.

The most common type of defect in the nematic phase is the disclination line defect, characterized by an integer or half-integer index defined by the number of turns of the director field associated with taking a circuit about the line. The phenomenological approach to the investigation of line defects was described by Schopohl and Sluckin schopohl.n:1987.a ; schopohl.n:1988.a . They worked with the full order-parameter tensor $`Q_{\alpha \beta }`$ in the framework of Landau-de Gennes theory, and applied this to the structure of the $`\pm \frac{1}{2}`$ disclination core. The *full* order parameter $`Q_{\alpha \beta }`$ was used a) to avoid *divergent* terms in the elastic energy; b) to take into account possible *biaxiality* of the defect core. Then, the general idea was implemented for the particular problem of the core structure of the $`\pm 1`$ strength defect: Biscari and Virga biscari.p:1997.a  and Mottram and Hogan mottram.nj:1997.a  solved the equations for the order tensor and obtained order parameter and biaxiality profiles. They used truncated expansions of the free energy density, which helped them to obtain analytical solutions to the problem. Finally, Sigillo et al sigillo.i:1998.a  considered the $`+1`$ disclination problem in the spirit of Maier-Saupe mean-field theory.

In this paper we use the approach proposed by Schopohl and Sluckin to obtain the order parameter profiles near the disclination line. We also recapitulate the results of the other theories in order to discuss the advantages and disadvantages of the proposed models. We compare these predictions with the results of computer simulation using a simple molecular model. Computer simulation is a well-established method of relating bulk elastic coefficients 1988.a ; 1990.a ; 1992.a ; 1996.g ; stelzer.j:1995.a ; stelzer.j:1995.b ; stelzer.j:1997.a  and more recently surface anchoring strengths stelzer.j:1997.b ; 1999.b , to molecular structure and interaction parameters. An early study of disclination line defects hudson.sd:1993.a  involved Monte Carlo simulations of rod-like molecules (hard spherocylinders) confined in cylindrical geometry with boundary conditions chosen to stabilize the chosen director field far from the defect. This work concentrated on the disclination line defect of strength $`\frac{1}{2}`$, in cylinders of radius 2–3 times the molecular length. For short rods, smooth variation of the order parameters with position was observed, and the defect core retained axial symmetry. For longer rods, having much larger bend elastic constants, this symmetry was broken, and various microscopically ‘escaped’ structures were seen. There was also evidence of metastability, and non-convergence of structures from different starting configurations, which the authors attempted to resolve by free energy calculations. More recently, a study of two-dimensional models has been carried out Dzubiella:1999.a . These studies support the view that the disclination core is one or two molecular lengths across.

The purpose of this paper is to present simulation results for a model of hard particles in a cylindrical pore, following closely the approach of Ref. hudson.sd:1993.a , and compare with the theories just mentioned. We study the simplest example of biaxial molecular arrangement, namely the disclination defect of strength $`+1`$ corresponding to a uniform, cylindrically symmetric, splay deformation of the director far from the core. This is imposed by choosing homeotropic anchoring conditions at the cylinder pore walls.

We take care to check the equilibration of our simulations, and note that quite extensive runs are required to ensure this. We also pay particular attention to studying the effects of pore size, with radius varying from 2 to 5 times the molecular length, so as to eliminate the effects of the walls on the defect core. More precisely, our aim is to restrict these effects to those which arise from the planar radial far-field structure, rather than from density and order parameter variations in the immediate vicinity of the walls. It is well-known that the pore radius plays a critical role in determining the stability of different nematic structures in this geometry. An analysis of the elastic free energy shows that the ‘escaped radial’ structure, in which the director bends over to become parallel with the symmetry axis in the core region, is the stable structure for large pore radius cladis.pe:1972.a ; kleman.m:1983.a . In analyzing experiments on nematic liquid crystals in cylindrical pores, Crawford et al. crawford.gp:1991.a ; crawford.gp:1992.a  show that at small pore radii, the escaped structure (which will also, in general, have point defects along the disclination line) is not the most stable; they discuss the planar radial and planar polar configurations, finding the latter to be stable for the parameter values that they survey.

A unified presentation of the different theories is given in Sec. II; the model, and simulation techniques, are set out in Sec. III. The observed structures are described in Sec. IV, along with theoretical fits to the order parameter profiles. A discussion of the results, and the validity of the theories, together with some concluding remarks, appear in Sec. V.

## II Phenomenological models

In the framework of the continuum theory, the system can be described by the Landau-de Gennes free energy density degennes.pg:1995.a :

$$F=\kappa \left|𝐐\right|^2+\sigma \left(𝐐\right),$$
(1)
where $`\sigma \left(𝐐\right)`$ is a function of the invariants of $`𝐐`$, the symmetric order tensor. This is the suitably-normalized traceless part of the tensor $`𝐌`$ of the second moments of the molecular orientational distribution function $`f`$

$`𝐌`$ $`=`$ $`{\displaystyle f\left(𝐮\right)𝐮𝐮d\mathrm{\Omega }},`$ (2)
$`𝐐`$ $`=`$ $`{\displaystyle \frac{3}{2}}𝐌{\displaystyle \frac{1}{2}}𝐈,`$
where $`𝐈`$ is the unit tensor and the integration is over the unit sphere.

We exclude the escape of the director in $`z`$ direction, and our boundary conditions also prevent spiral configurations, so the eigenvectors of $`𝐐`$ coincide with the unit vectors of the cylindrical geometry, namely $`𝐞_\rho `$, $`𝐞_\theta `$ and $`𝐞_z`$, respectively, in the radial, tangential, and axial directions. Since the eigenvalues of $`𝐐`$ are not independent, the tensor may be expressed in terms of two independent parameters: the order parameter $`S`$ and biaxiality $`\alpha `$ biscari.p:1997.a :

$`𝐐`$ $``$ $`Q_{\rho \rho }𝐞_\rho 𝐞_\rho +Q_{\theta \theta }𝐞_\theta 𝐞_\theta +Q_{zz}𝐞_z𝐞_z`$ (3)
$`=`$ $`S𝐞_\rho 𝐞_\rho +\left({\displaystyle \frac{1}{2}}S+{\displaystyle \frac{3}{2}}\alpha \right)𝐞_\theta 𝐞_\theta `$
$`+`$ $`\left({\displaystyle \frac{1}{2}}S{\displaystyle \frac{3}{2}}\alpha \right)𝐞_z𝐞_z.`$
Minimizing the full free energy functional (1) with respect to the $`S,\alpha `$ and taking into account that $`S=S\left(\rho \right)`$, $`\alpha =\alpha \left(\rho \right)`$, where $`\rho `$ is the distance from the cylinder axis, we obtain the following Euler-Lagrange equations:

$`\left(\rho S^{}\right)^{}3\rho ^1\left(S\alpha \right)3g_1(S,\alpha )`$ $`=`$ $`0,`$ (4)
$`\left(\rho \alpha ^{}\right)^{}+\rho ^1\left(S\alpha \right)g_2(S,\alpha )`$ $`=`$ $`0.`$
Here $`g_1(S,\alpha )=\frac{2}{9}\kappa ^1\sigma /S`$, $`g_2(S,\alpha )=\frac{2}{9}\kappa ^1\sigma /\alpha `$, and the prime denotes differentiation with respect to $`\rho `$.

The system of equations (4) should be solved with corresponding boundary conditions. The boundary conditions at $`\rho =0`$ can be obtained if we look for the solutions close to the centre of the disclination line. Indeed, in the limit $`\rho 0`$, we can neglect the functions $`g_1,g_2`$ in equations (4). Then, seeking the solution as an expansion in powers of $`\rho `$, one can obtain that in the region close to the centre of the disclination line the solutions are

$$S=S_0+3\gamma \rho ^2,\alpha =S_0\gamma \rho ^2,$$
(5)
where $`S_0`$ and $`\gamma `$ are constants. Therefore, at $`\rho =0`$, we have

$$S^{^{}}|_{\rho =0}=0,\alpha ^{^{}}|_{\rho =0}=0.$$
(6)
We also assume that, far away from the disclination core, we have uniaxial nematic, and the surface at $`\rho =R`$ provides the order parameter $`S_s`$:

$$S|_{\rho =R}=S_s,\alpha |_{\rho =R}=0.$$
(7)
Some general properties of the equations (4), regardless of the explicit form of the functions $`g_1`$, $`g_2`$, can help us to fit the simulation results. The simulations provide the value of the order parameter on the disclination line $`S_0`$ which we can also derive analytically. Indeed, the solutions (5) imply that $`S+3\alpha |_{\rho 0}=\text{const}`$. The condition that the energy is bounded requires also $`S_0=\alpha _0`$. Therefore, from (4) we obtain the implicit equation for the order parameter value on the disclination line:

$$g_1(S_0,S_0)+g_2(S_0,S_0)=0.$$
(8)
Now we consider the particular models.

### II.1 Full free energy expansion

A widely used formula for $`\sigma `$ is Landau-de Gennes’:

$$\sigma _{LG}=a\mathrm{Tr}𝐐^2b\mathrm{Tr}𝐐^3+c\left[\mathrm{Tr}𝐐^2\right]^2,$$
(9)
where $`a`$ is assumed to depend linearly on the temperature, whereas positive constants $`b`$, $`c`$ are considered temperature independent. For this free energy, the *uniaxial nematic* state is stable when $`b^2>24ca`$ with the degree of orientational ordering

$$S_b=\frac{b}{8c}\left(1+\sqrt{1\frac{64ca}{3b^2}}\right).$$
(10)
Taking into account the parametrization (3) the suitably rescaled potential (9) may therefore be rewritten as

$`\sigma _{LG}=k\{S^2({\displaystyle \frac{1}{2}}S^2{\displaystyle \frac{2}{3}}S(S_b+S_u)+S_bS_u)`$
$`+\alpha ^2[3S_uS_b+6(S_u+S_b)S+3S^2+{\displaystyle \frac{9}{2}}\alpha ^2]\},`$ (11)
so that the turning points for the uniaxial nematic phase with $`\alpha =0`$ occur at $`S=S_u,S_b`$.

The functions $`g_1,g_2`$ then read

$`\lambda ^2g_1`$ $`=`$ $`S_uS_bS+\left(S_u+S_b\right)\left\{3\alpha ^2S^2\right\}+S\left\{3\alpha ^2+S^2\right\},`$
$`\lambda ^2g_2`$ $`=`$ $`3S_uS_b\alpha +6\left(S_u+S_b\right)S\alpha +3\alpha \left\{3\alpha ^2+S^2\right\},`$
where $`\lambda ^2=\frac{9}{4}\kappa /k`$ is the characteristic length.

The equation (8) provides us with the values of the order parameter and biaxiality on the disclination line:

$$S_0=\alpha _0=\frac{1}{2}S_u.$$
(12)
### II.2 Biscari and Virga approach

In order to obtain analytical expressions for the order parameter profile, Biscari and Virga biscari.p:1997.a  used a quadratic approximation to the full free energy $`\sigma `$, expressed in terms of $`S`$ and $`\alpha `$:

$$\sigma _{BV}=k\left\{\eta \left(SS_b\right)^2+\alpha ^2\right\}$$
(13)
with the following expressions for the functions $`g_1,g_2`$:

$$g_1=\frac{\eta }{\lambda ^2}\left(SS_b\right),g_2=\frac{1}{\lambda ^2}\alpha .$$

The order parameter at the center of the disclination line then reads

$$S_0=\frac{\eta }{1+\eta }S_b.$$
(14)
### II.3 Mottram and Hogan approach

Another model, considered by Mottram and Hogan mottram.nj:1997.a  uses a quartic potential in $`S`$ but retains a quadratic potential in $`\alpha `$

$$\sigma _{MH}=k\left\{\eta S^2\left(\frac{1}{2}S^2\frac{2}{3}S\left(S_b+S_u\right)+S_bS_u\right)+\alpha ^2\right\}$$
(15)
This potential provides the following functions $`g_1,g_2`$:

$$g_1=\frac{\eta }{\lambda ^2}S\left(SS_u\right)\left(SS_b\right),g_2=\frac{1}{\lambda ^2}\alpha ,$$

and the following value of the order parameter on the disclination line:

$$S_0=\frac{1}{2}\left\{S_u+S_b\pm \sqrt{\left(S_uS_b\right)^24\eta ^1}\right\}$$
(16)
### II.4 Maier-Saupe approach

Another mean-field approach to the description of the disclination line was considered by Sigillo et al sigillo.i:1998.a . They applied Maier-Saupe theory and obtained expressions for the functions $`g_1,g_2`$. The Maier-Saupe theory has an intermolecular potential strength $`U`$ as the only adjustable parameter, determining the correct value of the bulk order parameter $`S_b\left(U\right)`$. The value of the order parameter $`S_0`$ on the disclination line is then fixed and determined by the value of $`S_b`$.

## III Molecular model and simulation methods

The molecules in this study are modelled as hard ellipsoids of revolution of elongation $`e=A/B=15`$, where $`A`$ is the length of the major axis and $`B`$ the length of the two equal minor axes. The phase diagram and properties of this family of models are well studied frenkel.d:1984.a ; frenkel.d:1985.a ; 1993.g ; 1993.j ; 1992.a . Units of length are chosen such that $`AB^2=1`$, making the molecular volume equal to that of a sphere of unit diameter. It is useful to express the density as a fraction of the close-packed density for perfectly aligned hard ellipsoids, assuming an affinely-transformed face-centred cubic or hexagonal close packed lattice; in reduced units $`\varrho _{\mathrm{cp}}^{}=\varrho _{\mathrm{cp}}AB^2=\sqrt{2}`$. Henceforth the asterisks denoting reduced quantities will be omitted.

The molecules were confined within a cylinder of radius $`R`$, and height $`H`$, with periodic boundary conditions applied in the $`z`$ direction (the symmetry axis). The cylindrical confining walls are defined by the condition that they cannot be penetrated by the *centres* of the ellipsoidal molecules. Packing considerations generate *homeotropic* ordering at the surface, i.e. the molecules prefer to orient normal to the interface, without the need to apply an explicit ordering field. The properties of such surfaces in planar geometry were investigated previously 1999.b . For these particles, the isotropic-nematic phase transition occurs at quite low density, $`\varrho /\varrho _{\mathrm{cp}}0.2`$. Temperature is not a significant thermodynamic quantity in this model; we set $`k_\mathrm{B}T=1`$ throughout.

Monte Carlo simulations were carried out for systems of the following cylinder radii: $`R/A=2.08,2.67,4.00,5.33`$. The corresponding numbers of ellipsoids were $`N=3500,6000,13000,22000`$. The height $`H/A`$ was in the range 2.6–2.8 in all cases. All the simulations were conducted at the same state point used in the earlier study 1999.b , corresponding to a bulk pressure $`P=2.0`$ in the above reduced units; this is significantly higher than the isotropic-nematic coexistence pressure $`P1.49`$ 1996.f . No external fields were applied, and conventional Monte Carlo moves were employed, with translational and rotational displacements chosen to give a reasonable acceptance rate.

The simulation results were analyzed to give a density profile, averaged over cylindrical shells of width $`0.125B`$$`0.25B`$, and an order tensor profile obtained by averaging

$$Q_{\alpha \beta }(\rho _i)=\frac{1}{n_i}\underset{j=1}{\overset{n_i}{}}\left\{\frac{3}{2}u_{j\alpha }u_{j\beta }\frac{1}{2}\delta _{\alpha \beta }\right\}\alpha ,\beta =\rho ,\theta ,z$$

where there are $`n_i`$ molecules present in bin $`i`$; $`\delta _{\alpha \beta }`$ is the Kronecker delta. The axis system is resolved as before into radial ($`\rho `$), tangential ($`\theta `$), and axial ($`z`$) components, for the purposes of accumulating these functions. Therefore, the tensor components $`Q_{\alpha \beta }`$ have been averaged over global rotations of the system (positions and orientations together) about the symmetry axis, as well as global translations in the $`z`$ direction. Diagonalizing this tensor, for each bin, gives three eigenvalues and three corresponding eigenvectors. This procedure allows us to test for (a) the planar radial structure, for which the eigenvectors lie along the $`\rho ,\theta ,z`$ coordinate axes, independent of $`\rho `$, and the eigenvalues are $`Q_{\rho \rho }(\rho )`$, $`Q_{\theta \theta }(\rho )`$, $`Q_{zz}(\rho )`$); and (b) the escaped radial structure, for which one of the eigenvectors points along the $`\theta `$ direction for all values of $`\rho `$, and the other two lie in the $`\rho ,z`$ plane, changing their orientation as $`\rho `$ varies.

The above tensor is not sensitive to any breaking of axial symmetry: instead, it gives a cylindrical average. To give a full representation of the positional variation of the orientational order, we must retain the full spatial dependence of $`𝐐`$. We find it convenient to calculate

$$S_{\alpha \beta }(𝐫)=\frac{1}{n_{}}\underset{j}{}u_{j\alpha }u_{j\beta }\alpha ,\beta =x,y,z$$
(17)
where the sum is conducted over the $`n_{}`$ molecules found to be in a neighbourhood $``$ of the chosen point $`𝐫`$. The eigenvectors of this $`𝐒`$ tensor are the same as those of the corresponding $`𝐐`$; the eigenvalues are linearly related and are non-negative, so the tensor may be visually represented as a spheroid, whose principal axis directions and lengths are given by the eigenvectors and corresponding eigenvalues 1989.a . Then, a well-ordered uniaxial region appears as an elongated, prolate, shape; the uniaxial core of a disclination defect of strength $`+1`$ would appear as an oblate spheroid with its symmetry axis along the axis of the cylinder; disordered regions would correspond to a spherical shape, and so on.

## IV Simulation results

Most results were obtained from starting configurations in which the ellipsoids were positioned randomly, avoiding overlaps, but perfectly aligned in the planar radial configuration. Typical equilibration times for these systems were around $`7\times 10^5`$ MC sweeps (one sweep is one attempted move per particle), and this was monitored through the orientation profiles defined above. The development of a biaxial core in the system with $`R/A=2.67`$ is shown in Fig. 1; the other systems evolved in a similar way. Following equilibration, production runs of approximately $`5\times 10^5`$ MC sweeps were undertaken.

Checks for convergence were carried out, starting from configurations in which molecular orientations as well as positions were disordered. The time evolution of the $`R/A=2.67`$ system is shown in Fig. 2. Note that the ordering remains uniaxial at the boundary, and propagates in towards the centre; as the ordering reaches the centre, the biaxial disclination core develops. An equilibration period of $`10^6`$ sweeps was needed to reach the equilibrium planar radial structure, and a subsequent production run of $`2\times 10^5`$ sweeps yielded identical results to those obtained from the planar radial starting point. For the larger cylinder radii, similar behaviour was observed, but the timescale for propagation of the order from the boundary to the centre was correspondingly longer. In each case, a run was conducted which was of sufficient length to confirm that the correct structure was becoming established.

In no case was the escaped radial structure formed during these runs. We have carried out some preliminary tests, for the pores of smaller radius, $`R/A=2.08,2.67`$, in which the escaped radial structure was stabilized by the application of a uniform orienting field favouring alignment along the cylinder axis. Following removal of this field, the planar structure was seen to be recovered on a simulation timescale of $`5\times 10^5`$ sweeps. Thus, the planar structure is thermodynamically stable for these cases, and it appears to be at least metastable on the timescales of our simulations for the larger pores.

In Fig. 3 the order parameter profiles are presented for all the different cylinder sizes. The biaxiality profile $`\alpha (\rho )=\frac{1}{3}(Q_{\theta \theta }Q_{zz})`$ indicates the extent of the core region. For the smallest two cylinders, it is clear that the walls act to deform the disclination core. For $`R/A=2.08`$, the core region is enlarged, while for $`R/A=2.67`$ it is compressed. The results for the two largest cylinders are almost indistinguishable.

In Fig. 4 we show the order tensor variation in the $`xy`$ plane (no variation with $`z`$-coordinate was observed) for the core region in the $`R/A=4.00`$ case. The $`𝐒`$-tensor spheroid of eqn (17) is plotted at a number of points $`(x_s,y_s)`$ which lie on circles centred on the cylinder axis, separated by $`0.25A`$, with the neighborhood $``$ of each point defined to include all molecules whose centres $`(x,y,z)`$ satisfy $`\sqrt{(xx_s)^2+(yy_s)^2}<0.25A`$. The tensors are averaged over a run of $`10^5`$ sweeps. The figure indicates that the core is actually best described as a pair of disclinations each of strength $`+1/2`$, symmetrically arranged at $`\rho /A`$1–1.25 from the cylinder axis. Very similar results are seen for all the three largest pores, $`R/A=2.67,4.00,5.33`$, while for the smallest pore, $`R/A=2.08`$, the defects are slightly further out, at $`\rho /A1.5`$. There is a small region of almost unperturbed uniaxial nematic liquid crystal around $`\rho 0`$, with the director perpendicular to the cylinder axis. We note that similar defect pairs are seen in the two-dimensional simulations reported in Ref. Dzubiella:1999.a .

The order tensor profiles reported in Fig. 3 are properly regarded as axial averages of structures which are not themselves axially symmetric on the timescales of the simulation. This is why the lowest eigenvalue in the profiles of Fig. 3 adopts similar values inside and outside the core: it corresponds to the eigenvector along the cylinder axis, and is unaffected by the axial averaging. We discuss in section V the validity of carrying out such an axial average, simply noting here that it is the most straightforward way of comparing our simulation results with theories which assume cylindrical symmetry.

We have fitted the results for $`R/A=4.00`$ using the phenomenological theories of section II. In each case, the boundary-value problem of eqns (4), (6), (7) , was solved using the *relaxation method* press.wh:1992.a . In Figs 5, 6, the fitting curves and the original simulation data are shown. Note, that to fit the simulation data we first determind the value of the parameter $`\eta `$ or $`S_u`$ using equations (12), (14), (16). Then we performed suitable rescaling, changing the factor $`\lambda `$.

## V Discussion

One can see in Figs 5, 6, that the slope of the fitting curves in general reflects the structure of the core: the center of the core is strongly biaxial, extending over a few units of length, as shown by the splitting of the eigenvalues $`Q_{\theta \theta }`$, $`Q_{zz}`$, and hence the non-zero biaxiality parameter $`\alpha (\rho )`$.

At the same time, the difference between the descriptions is also evident. Biscari and Virga’s theory gives an incorrect shape of $`S\left(\rho \right)`$ for small values of $`S`$ (and hence, small values of $`\rho `$ here). This is fairly predictable, since a quadratic expansion of the free energy density $`F\eta \left(SS_b\right)^2`$ was used, which is valid only for small deviations of the order parameter from the *bulk* value $`S_b`$.

The more sophisticated form of the free energy used by Mottram and Hogan (up to fourth order in $`S`$) corrects the slope of the curve for the order parameter $`S=Q_{\rho \rho }`$. However, the *biaxial* part, $`Q_{\theta \theta }Q_{zz}`$, still has only qualitative agreement with the simulation results, which is probably because of the *quadratic* approximation to the biaxial part of the free energy. This is particularly apparent in the $`\alpha (\rho )`$ profile, in which the biaxiality extends far beyond the core radius which would be deduced from $`Q_{\rho \rho }`$.

In spite of having only one adjustable parameter, the Maier-Saupe theory of Sigillo et al. gives a very realistic description of the disclination core structure. However, it predicts slightly lower core biaxiality than obtained in the simulations.

The most accurate fitting we obtained used the full expansion of the free energy in powers of the order parameter tensor (9), as shown in Fig. 6. This form, originally used by Schopohl and Sluckin, manages to reproduce the overall extent of the biaxial region quite well, while fitting the limiting behaviour of both $`\alpha (\rho )`$ and $`S(\rho )`$ as $`\rho 0`$, and the magnitude of $`S`$ in the bulk. Therefore, we can conclude, that all the terms in the free energy expansion (or, at least up to the fourth order in the order parameter tensor) should be taken into account if one want to make the correct *quantitative* description of the disclination core structure.

Fig 4 shows that the core structure seen on the simulation timescale is not, in fact, cylindrically symmetrical, and this raises questions regarding the validity of averaging the order tensor over rotations about the cylinder axis. Such averaging may arise in a natural way. We see a small amount of rotation of the defect pair around the axis, of the order of 10–15 degrees, as the core structure evolves in time, during the $`5\times 10^5`$-sweep simulation runs conducted here. If we roughly equate this to a real-time period of the order of nanoseconds, it seems possible that significant rotation will occur on the experimental timescale. The situation is further complicated when one considers correlations along the $`z`$-direction; the periodic box employed here is quite small, but use of a much longer cylindrical pore might reveal some twisting of the separate $`+1/2`$ defect lines about the cylinder axis. Both effects would result in the kind of cylindrical averaging which we have carried out, but we can currently only speculate on this.

A second question concerns the distance between the $`+1/2`$ defect lines. We do not observe separation of the defect pair, either as the simulations proceed in time, or as we study larger pores. Nonetheless, it is quite possible that the pore walls are confining the defect pair to the vicinity of the axis, and that they might separate if a much larger system were simulated. We note that, if they were to completely separate and approach the pore walls, the result would be the planar polar structure discussed elsewhere crawford.gp:1991.a ; crawford.gp:1992.a . It would be interesting to calculate free energies as a function of the defect separation to investigate this.

Finally, we may expect the escaped radial structure to become increasingly favoured as the pore radius increases; it is clearly not the most stable structure for the smallest pores studied here. Once more, free energy calculations are needed to make a proper comparison. Further work on these aspects is in progress.

## Acknowledgements

The advice of Nigel Mottram is gratefully acknowledged, and the referee’s comments on the first version of this paper were very helpful. This research was supported by EPSRC, and by the University of Bristol Department of Physics. DA also acknowledges support through the grants PSU082002 of the International Soros Science Education Program and ORS/99007015 of the Overseas Research Students Award. The authors of Ref. Dzubiella:1999.a  are thanked for providing a preprint.