A Strong Coupling Expansion

Strong Coupling Expansion of the Extended Hubbard Model with Spin-Orbit Coupling

Abstract

We study the strong coupling limit of the extended Hubbard model in two dimensions. The model consists of hopping, on-site interaction, nearest-neighbor interaction, spin-orbit coupling and Zeeman spin splitting. While the study of this model is motivated by a search for topological phases and in particular a topological superconductor, the methodology we develop may be useful for a variety of systems. We begin our treatment with a canonical transformation of the Hamiltonian to an effective model which is appropriate when the (quartic) interaction terms are larger than the (quadratic) kinetic and spin-orbit coupling terms. We proceed by analyzing the strong coupling model variationally. Since we are mostly interested in a superconducting phase we use a Gutzwiller projected BCS wavefunction as our variational state. To continue analytically we employ the Gutzwiller approximation and compare the calculated energy with Monte-Carlo calculations. Finally we determine the topology of the ground state and map out the topology phase diagram.

I Introduction

First introduced in 1963Hubbard (1963); Gutzwiller (1963); Kanamori (1963), the Hubbard model is regarded by many to be the simplest possible Hamiltonian which captures the essential physics of many-body systems with strong interactions. Owing to its simplicity, the Hubbard model has been used for decades to describe a variety of systems. Its applications have ranged from antiferromagnetismAnderson (1963) to the treatment of the metal-insulator transitionMott (1968); Hubbard (1964).

The standard Hubbard model contains hopping on the lattice sites (with hopping amplitude denoted by ) and on-site repulsion in the form of an energy penalty whenever two electrons are on the same site. It has been studied extensively in two dimensions on the square lattice mainly in the context of the high T cupratesScalapino (2006). The strongly interacting limit of this model at half filling is the Heisenberg Hamiltonian with antiferromagnetic coupling . Close to half filling (but not quite there) an appropriate approximation for the Hamiltonian is the t-J model where some hopping is allowed. This model is challenging to deal with since it contains both quadratic (kinetic) terms and quartic (spin) terms. Over the years many approaches where developed for the study of this strong interaction physics both analytically and numericallyLee et al. (1998). However, it is fair to say that at arbitrary doping (where is the hole density measured from half filling) any treatment uses some additional approximations.

In this paper we revisit the strong coupling limit of the Hubbard model with two additional terms (and hence the terminology ‘extended Hubbard’). The first is spin-orbit coupling which results in additional quadratic terms. These can be regarded as spin dependent/spin flip hopping processes. The second is off-site electron-electron interaction. We use this term to emulate the effect of the full Eliashberg method on the four fermion vertexOnari et al. (2006). In other words, instead of renormalizing the interaction vertex by the electron polarization bubble we add nearest neighbor attraction, . This effective attraction is appropriate whenever the polarization bubble of the fermions is maximal at Scalapino (1995).

The model we consider here is motivated by recent interest in topological superconductors and their promise to support Majorana fermions as bound excitations in vortex cores. The realization of topological superconductors has been proposed in different setups such as semiconductor heterostructuresFu and Kane (2008); Sau et al. (2010); Alicea (2010) and devices containing nanotubesLutchyn et al. (2010); Oreg et al. (2010); Cook and Franz (2011); Klinovaja et al. (2012). In these heterostructures the topological superconductivity is a result of spin-momentum locking provided by the spin-orbit coupling and tendency for pairing is induced through proximity to a superconductor. Inspired by the above advances we have proposed a model in which superconductivity arises from interactions (rather than proximity) in a system with spin-orbit couplingFarrell and Pereg-Barnea (2013a).

While the above model was studied in the weak coupling limitFarrell and Pereg-Barnea (2013a), superconductivity occurs at intermediate to strong coupling. In this Paper we lay the foundation of a strong coupling study of the extended Hubbard model in the presence of off-site interaction as well as spin-orbit coupling. We believe that this versatile model is useful beyond the scope of our specific application. For example, it might be useful in describing complex oxide heterostructure interfaces where the inversion symmetry is broken due to the superlattice. The broken inversion symmetry gives rise to Rashba spin-orbit coupling at the interface. Further, we feel that the methods developed here may be of use in other applications. One example of this might be the Kane-Mele-Hubbard model, which, among other things, is relevant to studies of how interactions effect topological band structuresRachel and Le Hur (2010); Reuther et al. (2012). It is therefore desirable to get a handle on the strong coupling limit of our model and this is the purpose of the present work.

To achieve the strong coupling limit of our model we assume that the interaction part of the Hamiltonian is large compared to the kinetic and spin-orbit coupling terms and develop a strong-coupling expansion. When this expansion is truncated at the second order and specialized to electron densities close to half filling it is a generalization of the model with Dzyaloshinskii-Moriya and compass anisotropy in the spin interaction. In this regime we employ the Gutzwiller approximationGutzwiller (1963); Zhang et al. (1988); Zhang (2003); Vollhardt (1984) to study the resulting model. This approximation has been successful in describing d-wave pairing and superconducting phase fluctuationsParamekanti et al. (2001, 2004) in the standard Hubbard model in two dimensions.

The Gutzwiller approximation consists of two stages. In the first a variational wavefunction is generated. When pairing is present, the appropriate variational wavefunction is a Gutzwiller projected BCS wavefunctionParamekanti et al. (2001); Sorella et al. (2002); Lugas et al. (2006). This wavefunction may contain any mean-field like orders (such as density waves, superconductivity or antiferromagnetism) and the Gutzwiller projection builds the strong interactions into it by eliminating any doubly occupied configurations. At this point one should evaluate the variational energy as the expectation value of the Hamiltonian with respect to the projected mean field wavefunction and minimize it with respect to the mean field order parameters.

The minimization of the variational energy in the many-body projected state is not trivial due to the projection operation and can not be done analytically. The evaluation should be done repeatedly until the energy is minimized and this can be done numerically using the Monte-Carlo techniqueFarrell and Pereg-Barnea (2013b). In the current work we choose to take a different approach and proceed analytically by making another approximation. Assuming that the charge (holes) is distributed uniformly on the lattice we may estimate the effect of the Gutzwiller projection on the various components of the variational energy. We may therefore renormalize the Hamiltonian parameters instead of performing the projection. This is the second stage of the Gutzwiller approximation which leads to a mean-field-like Hamiltonian whose parameters depend on the filling and can be analyzed in the standard way.

The rest of this article is organized as follows. In the next section we introduce the model and briefly present the generalized strong coupling expansion while leaving the details to Appendix A. Once we have obtained the effective Hamiltonian, we project it on a sub-space which is close to half-filling on the hole doped side and obtain a generalization of the model. We then study the strong coupling Hamiltonian first without spin-orbit coupling and a Zeeman term in section III, and then with finite spin-orbit coupling in section IV. We conclude with a mapping of the superconductor’s phase diagram according to the ground state topology.

Ii Strong Coupling Expansion

ii.1 The Model

The extended Hubbard model we consider is given by the following Hamiltonian on a two-dimensional square lattice:

(1)

where

(2)

is the tight binding kinetic energy and

(3)

where , (with a vector of Pauli matrices acting on the spin) and ( and are material parameters which describe the various spin-orbit coupling and Zeeman strengths). The three parameters in the spin-orbit coupling model above may originate from a variety of different sources. For example the parameters may come from a traditional spin-orbit model like the Rashba and Dresselhaus terms in Refs. [Sau et al., 2010,Alicea, 2010]. A second source may be parameters such as those used in the BHZ model, suitable for quantum wellsRothe et al. (2012); Lu et al. (2012); Guigou et al. (2011). Similarly, the Zeeman field parameter may be the result of a band gapBernevig et al. (2006), an external magnetic field or a magnetic field of a nearby ferromagnetic layerSau et al. (2010). As may come from a variety of sources we will ignore any orbital effects that could arise in some cases.

The interaction terms denoted by above are given by

(4)

where is the on-site repulsion strength and describes attraction between nearest neighbors.

For the purpose of making the strong coupling expansion we will need to express the Hamiltonian in real space. We therefore transform the spin orbit coupling to real space and divide it into two types of terms. , the Zeeman Hamiltonian will contain all on-site terms, proportional to :

(5)

The other terms in amount to hopping processes which act non-trivially on the spin. We write them in real space and combine them with the hopping into a generalized hopping matrix such that:

(6)

The reader should note that this model has been proposed not with a specific physical system in mind but for the sake of versatility. The work to follow is meant to motivate a search for a specific material with these properties. One possible candidate that may be described by this model is copper intercalated BiSe which, although a 3D material, develops a 2D-like Fermi surface at certain dopingLahoud et al. (2013).

ii.2 The strong coupling Hamiltonian

The strong coupling expansion distinguishes between the high interaction energy scales and other, quadratic terms in the Hamiltonian. Written in real space the quadratic terms are either on-site (chemical potential, Zeeman) or hopping. The on-site terms do not change the potential energy while the hopping terms do. Ideally, one would like to diagonalize the Hamiltonian but due to its quartic terms this can not be done analytically. Instead of diagonalizing we set out to block diagonalize the Hamiltonian. Since we are interested in the strong coupling regime we would like to block-diagonalize the Hamiltonian such that the interaction energy is constant at each block and work at the lowest interaction energy block. In other words - there exist a unitary transformation which eliminates terms which change the interaction energy from the Hamiltonian. Using the properties of the desired transformation (unitarity and interaction energy conservation) we can formally write down the transformation. In order to find a closed form, however, we resort to a power expansion in a small parameter of the order of the ratio between quadratic part of the energy and the interaction part of it. We are then able to find the transformation and the transformed Hamiltonian up to a given order.

Figure 1: Schematic diagram of the a hopping process contributing to .

Let us sketch the procedure which resembles the treatment of ref’s. [MacDonald et al., 1988] of the original Hubbard model. We denote the unitary transformation that block diagonalizes the Hamiltonian by and the Hamiltonian terms excluding the (constant) interaction and Zeeman energy by . In order to perform a series expansion for both the unitary transformation and the transformed Hamiltonian we arrange the terms in by their effect on the interaction energy. For example, the operator includes all terms which do not change the interaction energy at all. We further define the operators as a collection of the hopping terms which change the number of doubly occupied sites by and the number of nearest neighbor pairs the electron sees from to . This amounts to a change of to the interaction energy. As an example of one of these processes we refer the reader to Fig. 1. In the figure we illustrate one of the hopping processes involved in the specific operator . An electron starts on lattice site and hops to a nearest neighbour. In the process of this hop we see a double occupancy is created and so . Further, the electron that we are hopping begins on lattice site which has 4 nearest neighbours and then hops (to the right in the figure) to a site with 3 nearest neighbours. Thus this process has and After a few more steps we arrive at the expansion and present it up to second order:

(7)
(8)

where the primed sum excludes terms in which .

The results presented in Eq. (7) are general for any . Moreover, since the denominator of the above expression is a combination of and it is sufficient to have only one of them large compared with the band width. In our system one may expect that the on-site repulsion be large while the nearest neighbor attraction is small. Nevertheless, our formalism allows us to treat both interaction terms exactly and we do so.

To proceed we focus on the parameter regime (and is actually smaller than such that it does not induce any charge clustering). For this parameter regime the ground state of a system on the hole-doped side of half-filling before we account for or will lie in the subspace of real-space configurations with no doubly occupied sitesMacDonald et al. (1988); van Dongen (1994). We consider a complete subspace of no doubly occupied sites even though the effects of may lead so some phase separation of electrons and holes. We feel this to be an appropriate simplification to make in order to make progress with the problem. At half filling the second part of the Hamiltonian that is second order in the operators can be rearranged as a spin Hamiltonian which has been studied elsewhereFarrell et al. (2013). At, or slightly hole doped away from, half-filling the dominant contribution from the term of the form will see an initial number of nearest-neighbours and hop to a site with . All other terms involved in the second order contribution are either impossible or occur with a probability proportional to some power of , which is very small near half-filling. Thus although other terms that are of order exist in the second order term in Eq. (7), the term with denominator is the statistically dominant term. A picture of this process can be formed by looking at Fig. 1. This process corresponds to which corresponds to the hop shown in Fig. 1 and then one of the two electrons on the site to the right of site hopping back. Close to, but not at half filling hopping without a change to the interaction energy is possible and included in . This is the generalization of the well known model to our case and the subject of this paper. The model is given by:

(9)

where the spin vector is defined by and the spin coupling is given by

(10)

where is the nearest neighbor vector and is for and for . From the above matrix one can read off the Heisenberg, Dzyaloshinskii-Moriya and compass anisotropy terms:

(11)

Spin models similar to the one above have been employed in the study of spin-orbit effects in ultra-cold atoms and other systemsBanerjee et al. (2013); Radić et al. (2012); Cole et al. (2012). The focus of the remainder of this paper will be on studying a renormalized version of the above model.

Iii A Gutzwiller Approximation at Strong Coupling

In order to find the ground state of the strong coupling limit of our model we adapt the Gutzwiller approximationZhang et al. (1988); Zhang (2003); Vollhardt (1984) to the extended Hubbard model. The rational of this approximation is as follows. The Gutzwiller approximation is a variational method that uses a projected BCS mean-field wavefunction.

A projected BCS wavefunction has proved useful in the context of strongly correlated electron systems. It is obtained through projecting out all doubly occupied configurations from the BCS wavefunction. The resulting wavefunction, encodes both the tendency for developing order such as Cooper pairing and density waves due to its mean field starting point as well as the strong on-site repulsion of the model due to the Gutzwiller projection. The reader should note that more sophisticated generalizations of this wave function have been used in Monte Carlo calculations in the pastCapello et al. (2005, 2006)

The variational energy should be calculated and minimized with respect to the order parameters built into . However, the application of the Gutzwiller projection on the wavefunction makes the evaluation very complicated and not attainable analytically. One may then resort to numerical methods such as Monte Carlo integration in order to find Paramekanti et al. (2001); Sorella et al. (2002); Lugas et al. (2006). Monte Carlo methods have been refined in recent years to the point that they can be regarded as variationally exact. We take this numerical approach elsewhereFarrell and Pereg-Barnea (2013b). Here we employ an analytical approach, namely, the Gutzwiller approximation. The Gutzwiller approximation which was developed several decades agoZhang et al. (1988); Zhang (2003); Vollhardt (1984) is a method to approximately apply the Gutzwiller projection operator. First, one should note that the projection operator can be applied either on the wavefunction or on both sides of the Hamiltonian. This means that one may derive an effective Hamiltonian that acts on the non-projected mean-field wavefunction. The resulting energy is therefore an approximation of the expectation value of the original Hamiltonian with respect to the projected wavefunction. In this section we apply the Gutzwiller approximation to the extended Hubbard model, starting without spin-orbit coupling and then adding it in the next section.

Without spin orbit coupling and Zeeman field Eq. (9) reduces to:

(12)

where the only difference from the usual model is in the kinetic term. In this term we only include hopping processes which do not change the number of doubly occupied sites as well as the number of nearest neighbor pairs. We assume hole doping, in the vicinity of half filling such that and . We further assume that the system is completely unpolarized, i.e. . While this may be expected for a spin conserving Hamiltonian, it may not be locally accurate for a system with spin-orbit coupling and certainly will not be favorable when a Zeeman term is included. The variational energy is given by:

(13)

Note that since our projection will be applied to the Hamiltonian rather than the wavefunction, we do not need to specify the mean-field wavefunction just yet.

First let us estimate the change in the the expectation value of the hopping term due to the Gutzwiller projection and the definition of which is restricted to hopping processes which do not change the interaction energy. The kinetic part of the variational energy is therefore

(14)

where the operator projects out all possible spin orientations of the nearest neighbors of site except the one specified by , where takes a value of 1 (a value of zero) for site occupied (unoccupied) by an electron with spin . A formal definition of this operator is left to Appendix A. Let us consider the action of on any configuration that may appear in the BCS state. First removes all possible double occupancies. For a system close to half filling this projection means that the only configurations left in are ones with singly occupied or unoccupied sites. then moves an electron from an occupied site with nearest neighbors to an empty site with nearest neighbors. Let us compare this to the normal unprojected operator, , acting on just the BCS ground state. When this term acts, it finds a lattice site which is occupied by an electron of some spin and then moves it to a neighbor site that is not occupied by this particular spin.

We can assign relative probabilities for each of these two events. For the latter case the probability that the process occurs is just the probability that the lattice site the electron moves to is void of the particular spin we are moving. Assuming that the probability that any site is occupied is uniform then the probability of the particular move we just described is . For the case of the projected expectation value there are several concurrent requirements we must meet. First, we require the lattice site where the electron is moved to be empty. This gives a factor of . Second, the site we move the electron from must have exactly occupied nearest neighbors sites. The site we are moving the electron to is a nearest neighbor site and must be empty, therefore . The probability of having occupied nearest neighbor sites is then the binomial probability of having exactly successes in 3 trials with success probability , in other words . Finally, the site we move to must also have exactly occupied nearest neighbors, this again occurs with probability . Combining these three observations we find that the probability of the entire series of events is .

This means that

(15)

and we therefore redefine

(16)

We now turn to the quartic terms. Although our goal is to obtain something that is written as a spin-spin interaction we consider the double hopping processes that give rise to them as this language is more suitable for the projection. The interaction energy can be written as the expectation value of the quartic terms in and reads

(17)

where the primed sum is as defined before.

First we repeat that in systems with hole doping and after is applied to the wavefunction, there are only singly occupied or empty sites left in each configuration in the mean field state. This means that can never be , as it is impossible to reduce the number of doubly occupied sites. Therefore is either or . We will ignore the terms above, these terms require the movement of two electrons to two empty sites and so are a factor of less probable than their counterparts. Let us then consider only . This term creates a doubly occupied site and then destroys it. This can happen in two ways: (1) an electron hops to an occupied site and then one of these two electrons hops back to the original site or, (2) an electron hops to an occupied nearest neighbor and then one of the two electrons on this site hops to an empty next-nearest neighbor site. Keeping with our motive of simplicity, we will ignore the latter “three site” move as it is less probable near half filling. Taking all of these concerns into account we have

(18)

We now estimate the Gutzwiller renormalization for . First consider a non-projected process. In this case an electron hops to an occupied site and then one of the two electrons at the site hops back. In order for this to occur both sites must only be occupied by one spin. The probability of having a site only occupied by (say) an up spin next to a site which is only occupied by a down spin is . Now we consider the probability with all of the appropriate projections. The net result of is that an electron with an up spin hops from a site with nearest neighbors to a site with which is occupied by a down spin with nearest neighbors and then one of the two electrons now occupying this nearest neighbour site hops back to the original site. can range from to while can take values from to . The probability of this process is . We therefore approximate the interaction energy by

(19)

where we have omitted the quadratic terms which originate from rearranging the fermion operators.

The next step in the Gutzwiller approximation is to evaluate the variational energy (this time with a non-projected BCS wavefunction) and minimize it with respect to the order parameters. As this step is similar to the standard variational mean field procedure we de not give the details here. Some sample results and their comparison to a numerical Monte-Carlo evaluation are given in Appendix B.

Iv Finite Spin-Orbit Coupling

In this section we keep the parameters and finite. Keeping either or finite leads to spin-polarization of the electrons in the system. This inequality in the number of spin-up and spin-down electrons complicates the nature of the Gutzwiller approximation as the tacit assumption we made in the last section that is now violated. We therefore must estimate the Gutzwiller renormalization for processes involving different spin flavors separately.

We begin by considering the hopping term in Eq.(7) and estimating the effect of the Gutzwiller projection. The first projection we deal with is the restriction on the number of nearest neighbours permitted before and after an electron hops from one site to the other. As this part of the projection depends only on the site occupation and not the spin it results in a factor similar to the one we had before:

(20)

Next we treat the different parts of separately. Recall that we defined the generalized hopping matrix, , in Eq. (II.1) as having two parts. One includes simple hopping (which may have different amplitudes for different spins) and the other includes processes in which the spin is flipped during the hopping. We estimate the effect of the projection on each one of the processes by considering all possible ways this process can occur in some state and all possible ways they can occur in a Gutzwiller projected state Zhang et al. (1988); Ko et al. (2007); Chou and Lee (2012). This leads to the renormalized hopping matrix.

(21)

We can see that is the relative probability of a hop of an electron with spin . In the projected term we require an empty site in order to make the hop while in the unprojected term we simply need the destination site to not be occupied by a spin . Similarly we can understand . The projected process requires an empty site and hence the factor of in the numerator. On the other hand, the unprojected process requires that the site the electron hops to has no electron with the opposite spin. This accounts for a factor , however we must equally weight this with the probability for the process to happen in reverse, this is where the factor of and the square root come from.

To obtain the Gutzwiller factor of the quartic term in our strong coupling model at finite spin-orbit coupling we (again) neglect three site hopping and processes where the (in the operators) is zero. This leaves us with only terms that begin and end on the same site which can be written as: . Let us rewrite the coupling matrix as:

(22)

where we have defined the three couplings , , and and the function is 1 if and -1 if . Here we have already approximately taken the restrictions on the number of nearest neighbors permitted by a given hop by replacing the denominator by the Gutzwiller renormalization which has been defined in Eq. (III). We have not yet taken the occupancy requirements into account. We now consider the factors which are spin dependent to obtain individual renormalization factors for different matrix elements.

To carry out this renormalization one can group the spin Hamiltonian into four separate terms, through . Each one of these four terms describes a different physical process. The first term, , flips the spins on nearest neighbour sites (e.g. an up spin on site and a down spin on site are flipped to a down and an up spin respectively). The second term looks at the component of the spins on neighbouring lattice sites and takes the form . The third term looks for nearest neighbour sites with two up or two down spins and flips both spins. Finally, the fourth flips a spin on site while measuring the component of a spin on site . By considering the amplitude for each of these four physical processes in the projected and unprojected state, we can form Gutzwiller renormalization factors for each of them. We abdicate the details of this process to Appendix C, the result of performing this weighting of amplitudes is as follows

(23)

After making the replacements for each of the four terms we can rewrite the spin-contribution as . In the first two diagonal entries the and terms are renormalized by . The Dzjaloshinskii-Moryia term is renormalized by and the diagonal term by . Finally, the term is renormalized by . Note that although formally we have obtained we have followed [Ko et al., 2007] and sent in our calculations to restore rotational symmetry that is broken by our Gutzwiller approximation. This replacement is not motivated by the specific model we have taken (where SU(2) is already broken), but instead is made to ensure that our Gutzwiller approximation hasn’t broken any additional symmetry. The role of the factors is to approximately take into account on statistical grounds. There is no reason why , which does not break rotational symmetry, should effect and (or ) differently. The replacement is made in order to respect this fact. Note that the roman numeral subscripts on the renormalization factors should not be thought of as renormalizing specific couplings , , etc. nor with the (italic) site index ; rather they correspond to the four physical processes that can result when acts on a real space configuration of electrons. For example, renormalizes the spin-flip processes that occur as a result of the term in entries and aboveFarrell (2013).

The Gutzwiller renormalized spin interaction matrix reads:

(24)

where , , and finally . With this renormalized exchange matrix we have the following fully projected Gutzwiller Hamiltonian

(25)

we now turn our focus to a mean field study of this Hamiltonian.

iv.1 Mean-field decoupling of the projected Hamiltonian

We are now ready to treat the approximately projected Hamiltonian in mean field. In order to do so we define an auxiliary Hamiltonian from which the variational wavefunction can be found.

(26)

where , and are all functions to be determined. Note that we must allow for four channels of pairing (including triplet pairing as well) due to the spin asymmetry which is built into the model.

We perform a mean field decoupling of the Gutzwiller projected Hamiltonian and compare the result to Equation (26). In doing so we find that the free parameters built into this auxiliary model have the following forms

(27)

where , we have defined the effective field and we have defined the following parameters

(28)

In the above there are twelve free parameters that must be fixed through ensuring that the theory is self-consistent. These self-consistency conditions are as follows

(29)

where to make the equations compact we have used and and while . These twelve equations must be solved simultaneously with the particle number constraints and . This amounts to a solving a fifteen by fifteen system of non-linear equations. For clarity let us explicitly mention the hierarchy of notation in moving from Eq. (27) to Eq. (28) and finally to Eq. (29). We have used regular symbols, e.g. , to denote the dependent values in the auxiliary Hamiltonian. The parameters with the hat symbol, e.g. , are constants that enter the definition of dependent functions used to define the auxiliary Hamiltonian. Finally, we have used the tilde symbol, e.g. , for self-consistent parameters that are heretofore unknown until we solve the equations in Eq. (29).

Note that triplet superconductivity arises from the self-consistency equations whenever spin-orbit or Zeeman terms are present. With finite spin-orbit coupling the conditions in Eq. (29) cannot simply be solved by setting the triplet superconductivity term to zero. This stems from the non-trivial coupling between and , and , and and and that is in place as a results of the coupling . This is seen by inspecting the last term in each of the definitions in Equation (28). To see why this coupling refuses solutions where some of the parameters are identically zero consider looking for a solution where . The last two equations in Eq. (28) then become

(30)

and the relevant equations we must solve are

(31)

It is impractical to find analytic expressions for averages such as (they depend on eigenvectors of a complicated four by four matrix) however we can find these numerically. Looking at numerics and using some intuition we argue that very roughly and we have said above that . So in general either or the sum on the right of the second equation in (31) is nonzero and therefore cannot be zero.

We have solved the equations developed in the theory above numerically. Where these parameters are non-zero the general trend we see in the solutions is , , and . This equality of parameters with and directionality is not surprising as our model has no preference for the or direction.

We first present our results for and . In this special case of the spin-orbit coupling parameters there is not preferred spin direction and we have , and