text:
# Hydrodynamic Instability of the Flux-antiflux Interface in Type-II Superconductors

## Abstract

The macroturbulence instability observed in fluxline systems during remagnetization of superconductors is explained. It is shown that when a region with flux is invaded by antiflux the interface can become unstable if there is a relative tangential flux motion. This condition occurs at the interface when the viscosity is anisotropic, e.g., due to flux guiding by twin boundaries in crystals. The phenomenon is similar to the instability of the tangential discontinuity in classical hydrodynamics. The obtained results are supported by magneto-optical observations of flux distribution on the surface of a YBCO single crystal with twins.

The interest in the physics of the vortex state in type-II superconductors increased significantly after the discovery of high-$`T_c`$ superconductivity (HTS). A main reason for the renewed attention is the observation of many novel nontrivial phenomena occurring in the vortex matter of the HTS materials. The perhaps most dramatic of these phenomena is the turbulence instability of the vortex-antivortex interface, which was observed in 1–2–3 systems using magneto-optical (MO) imaging . It consists in the following. When magnetic flux is trapped in the superconductor and a moderate field of reverse direction is subsequently applied, a boundary of zero flux density will separate regions containing flux and antiflux. In some temperature and field range this flux-antiflux distribution can display unstable behavior characterized by an irregular propagation of the boundary where finger-like patterns develop. This contrasts strongly the regular propagation of the flux front during virgin field penetration, when only one flux polarity is present in the sample. It is clear that the instability, often called macroturbulence in the literature, cannot be understood within the frameworks of the critical state model or conventional models for flux relaxation .

An attempt to explain this remarkable behavior of the flux-antiflux interface was made in , where the instability was attributed to a thermal wave generated by local heat release in the vortex-antivortex annihilation. Unfortunately, this mechanism can hardly be accepted. Indeed, the vortex energy consists of two terms, where one is the magnetic part related to the magnetic field of the vortices. The other part is stored within the vortex core and represents the condensation energy. The magnetic energy of a vortex-antivortex pair is dissipated as Joule heat as they are getting close to each other, but before the annihilation takes place. This energy, which equals $`\stackrel{}{J}\stackrel{}{E}`$ where $`\stackrel{}{J}`$ is the current density and $`\stackrel{}{E}`$ is the electric field generated by the vortex motion, is dissipates not only at the interface, but in the bulk of the sample. Hence, only the release of the core energy in the process of annihilation is concentrated near the interface. A simple calculation shows that the core energy is much smaller than the magnetic part and will cause only a negligible rise of the sample temperature , thereby ruling out that thermal effects are responsible for the instability.

In this Letter we present an explanation focussing on the experimental fact that the instability was reported only for $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_{7\delta }`$ and other 1–2–3 single crystals. These materials are characterized by the existence of two twin boundary systems oriented orthogonally to each other in the ab plane. The vortex motion in such crystals is governed by a pronounced guiding effect , and the viscosity coefficient for motion along the twin boundaries should be much smaller than across them. This anisotropy gives rise to vortex motion with a velocity component normal to the Lorentz force. The vortices and antivortices are then forced to move towards each other along the interface, where the tangential component of the velocity becomes discontinuous. It is well known that hydrodynamical flow under such conditions can be unstable , and we show here that a purely hydrodynamic theory of the vortex-antivortex system with anisotropic viscosity can explain the origin of the macroturbulence.

Consider an infinite superconducting plate of thickness $`2d`$ with the external magnetic field $`\stackrel{}{H}`$ oriented parallel to the sample surface along the $`z`$-axis. The $`x`$-axis is perpendicular to the plate and $`x=0`$ in the center. Let $`H`$ first increase and then be lowered through zero to a negative value. Then two kinds of vortices will exist in the sample; one with field direction along the positive $`z`$-axis (vortices) and one directed oppositely (antivortices). From the symmetry of the problem it is sufficient to consider only the region $`0<x<d`$, and Fig. 1 shows schematically the distributions $`N_1(x)`$ and $`N_2(x)`$ of vortex and antivortex densities.

The densities $`N_1`$ and $`N_2`$ satisfy the continuity equation,

$$\frac{N_\alpha }{t}+\mathrm{div}(N_\alpha \stackrel{}{V}_\alpha )=0,\alpha =1,2$$
(1)
where $`\stackrel{}{V}_\alpha `$ are the vortex and antivortex velocities. In the regime of anisotropic viscous flow they are related to the Lorentz driving force by,

$$\eta _{ik}N_\alpha V_{\alpha k}=F_{Li},\stackrel{}{F}_L=\frac{1}{c}\stackrel{}{B}\times \stackrel{}{J},$$
(2)
$`\stackrel{}{B}=N_\alpha \stackrel{}{\mathrm{\Phi }}_0,\stackrel{}{J}={\displaystyle \frac{c}{4\pi }}\times \stackrel{}{B}.`$

Here $`\eta _{ik}`$ is the symmetrical tensor of the anisotropic viscosity and $`\mathrm{\Phi }_0`$ is the magnetic flux quantum. Eq. (2) can be rewritten as,

$$V_{\alpha i}=\frac{\mathrm{\Gamma }\mathrm{\Phi }_0^2}{4\pi }\gamma _{ik}\frac{N_\alpha }{x_k},$$
(3)
where $`\gamma _{ik}`$ is the dimensionless tensor of the inverse viscosity, $`\eta _{ik}^1=\mathrm{\Gamma }\gamma _{ik}`$, and the coefficient $`\mathrm{\Gamma }`$ is chosen so that the principal values of $`\gamma _{ik}`$ are unity and $`\epsilon `$, which satisfies $`0<\epsilon <1`$. The case of a strong guiding effect corresponds to $`\epsilon 0`$. Note that $`\mathrm{\Gamma }`$ increases rapidly with the temperature due to thermal depinning of the vortices.

To solve the problem one must also formulate boundary conditions at the vortex-antivortex interface. Generally the position of the interface, $`x=x_0`$, depends on the $`y`$-coordinate and time $`t`$ as it moves with a velocity $`\stackrel{}{U}`$. The first condition is that the total flux of vortices and antivortices through the interface vanishes,

$$N_1(\stackrel{}{V}_1\stackrel{}{U})_n+N_2(\stackrel{}{V}_2\stackrel{}{U})_n=0.$$
(4)
Secondly, both vortex polarities annihilate at the interface with the rate proportional to the product of their densities,

$$N_1(\stackrel{}{V}_1\stackrel{}{U})_n=RN_1N_2.$$
(5)
The parameter $`R`$ can depend on the vortex densities and velocities but, for simplicity, we assume here that it is a phenomenological constant. A similar model for the annihilation process was used in Ref.. However, contrary to , we consider the parameter $`R`$ to be defined by the microscopic Meissner current of individual vortices, which is much greater than the macroscopic currents $`J`$. Therefore, the relative velocity of annihilating vortices and antivortices is much greater than the hydrodynamic velocity $`V_\alpha `$. As a result, the region where vortices and antivortices coexist and annihilation takes place is very narrow, and can be represented by the surface $`x=x_0(y,t)`$. Finally, we take the magnetic induction to be zero, i. e.,

$$N_1=N_2,$$
(6)
at $`x=x_0(y,t)`$. This condition is also a direct consequence of Eq. (4) if the vortex and antivortex densities on the interface were equal to zero initially.

As a first step, we determine what the model predicts for the unperturbed distribution profiles $`N_1`$ and $`N_2`$. In this case the vortex-antivortex interface is the plane $`x=x_0(t)`$ which moves with the velocity $`U=\mathrm{d}x_0(t)/\mathrm{d}t`$. Since a simple solution with constant $`U0`$ could not be found, we consider only the stationary profile with $`U=0`$. With $`N_1/t=N_2/t=0`$, the Eqs. (1), (4)–(6) give,

$$N_1(x)=N_2(d)\left(\frac{x_0+d/2rx}{d+d/2rx_0}\right)^{1/2},$$

$$N_2(d)=H_0/\mathrm{\Phi }_0,$$
(7)
$$N_2(x)=N_2(d)\left(\frac{x+d/2rx_0}{d+d/2rx_0}\right)^{1/2},r=\frac{4\pi Rd}{\mathrm{\Gamma }\mathrm{\Phi }_0^2\gamma _{xx}}.$$

In the following, we assume $`r1`$ since the rate of the annihilation is fast and the viscosity is not small. Therefore, at the interface the vortex densities are relatively small, $`N_1=N_2N_2(d)r^{1/2}`$ while the spatial derivatives of $`N_1`$ and $`N_2`$ at $`x=x_0`$ are large,

$$N_2^{}=N_1^{}=N_1r/d,N_2^{\prime \prime }=N_1^{\prime \prime }=N_1(r/d)^2.$$
(8)
Then, also the velocities $`V_{\alpha x,y}(x=x_0)`$ are relatively high, since they are proportional to $`r^{1/2}`$ near the interface.

If one cannot provide fixed values of the vortex and antivortex densities $`N_1(0)`$ and $`N_2(d)`$ at the boundaries, the discussed stationary profiles cannot be realized. In practice the external magnetic field defines the density of the antivortices only. Because of this asymmetry and the vortex-antivortex annihilation, the total number of the vortices will decrease with time and the interface moves towards the middle of the sample. To simplify the problem, we assume the interface velocity $`U`$ to be much smaller than the vortex velocity $`V_\alpha `$, i. e. we consider the problem in the quasi-stationary regime.

To investigate the stability of the interface with respect to small perturbations, it is suitable to introduce the following dimensionless variables,

$`n_\alpha `$ $`=`$ $`N_\alpha /N_\alpha (x_0),\tau =t/t_0,t_0={\displaystyle \frac{\mathrm{\Gamma }\mathrm{\Phi }_0^2\gamma _{xx}}{4\pi R^2N_\alpha (x_0)}},`$ (9)
$`\xi `$ $`=`$ $`x/L,\zeta =y/L,L={\displaystyle \frac{\mathrm{\Gamma }\mathrm{\Phi }_0^2\gamma _{xx}}{4\pi R}}=d/r.`$ (10)
Normalization using the time-dependent $`N_\alpha (x_0(t))`$ is allowed here since we assume that the instability develops much faster than noticeable changes occur in $`N_\alpha (x_0)`$. We seek for perturbations in the vortex and antivortex densities of the form

$$n_\alpha =n_\alpha ^{(0)}+f_\alpha \mathrm{exp}[\lambda \tau +\mathrm{i}k\zeta +p_\alpha (\xi \xi _0(\tau ))].$$
(11)
The linearized boundary conditions should be written on the perturbed interface,

$$\xi =\xi _0(\zeta ,\tau )=\xi _0(\tau )+\delta \xi \mathrm{exp}(\mathrm{i}k\zeta +\lambda \tau ),$$
(12)
with the normal unit vector

$$\stackrel{}{\nu }=(1,\mathrm{i}k\delta \xi (\zeta ,\tau )).$$
(13)
It follows directly from Eq. (6) that

$$\delta \xi =(f_1f_2)/2.$$
(14)
Equations (1) give the expressions for the parameters $`p_1`$ and $`p_2`$. Substituting them and Eqs. (11), (14) into Eqs. (4), (5), we obtain two linear algebraic homogeneous equations relating $`f_1`$ and $`f_2`$. Demanding the determinant to vanish and omitting terms of the order of $`(Ut_0/L)^21`$, one obtains the following dispersion equation for the increment $`\lambda `$ at different wave numbers $`k`$,

$$\lambda =\mathrm{\Omega }^2ϵ\kappa ^22\mathrm{i}s\kappa 1b.$$
(15)
Here

$$\kappa =\frac{k|\alpha u|}{2},\alpha =\frac{\gamma _{xy}}{\gamma _{xx}},u=U\frac{t_0}{L},$$
(16)
$`ϵ=4\epsilon /(\alpha u)^2,s=\mathrm{sign}(\alpha u),2b=n_1^{\prime \prime }n_2^{\prime \prime },`$

and $`\mathrm{\Omega }`$ is a root with $`\mathrm{Re}\mathrm{\Omega }>0`$ of the equation

$`\mathrm{\Omega }^4+3\mathrm{\Omega }^3+\mathrm{\Omega }^2(ϵ\kappa ^2+22b)`$ (17)
$`\mathrm{\Omega }(2ϵ\kappa ^2+\mathrm{i}s\kappa +4b)3s\mathrm{i}\kappa =0.`$ (18)
Shown in Fig. 2 is the dependence of the increment $`\mathrm{Re}\lambda `$ on the dimensionless wave number $`\kappa `$ for different $`ϵ`$. The parameter $`b`$ is set to unity as this follows from the quasi-stationarity condition $`u0`$. The curves demonstrate that a positive increment exists, i.e. the planar interface becomes unstable, when $`ϵ<ϵ_c=0.019`$. The instability occurs at not too small values of the anisotropy and velocity of the flux-antiflux interface. Furthermore, it is characterised by a temporal scale given by the largest $`\mathrm{Re}\lambda `$, $`\lambda _m`$, which occurs at finite $`\kappa =\kappa _m`$. For $`ϵ`$ sufficiently small, $`\kappa _m`$ is much greater than unity and one finds from Eq. (17)

$$\lambda _m=1/(4\sqrt{ϵ})2,\kappa _m=1/\sqrt{2}ϵ^{3/4}.$$
(19)
Thus, we have identified that an instability occurs if $`ϵ<ϵ_c`$, which in dimensional notation is expressed as

$$\epsilon <ϵ_c\left[\frac{U\mathrm{tan}\theta }{2RN_1(x_0)}\right]^2$$
(20)
where $`\theta `$ is the angle between the direction of flux guiding and the flux-antiflux interface. The above analysis allows us to understand why the instability of the vortex-antivortex interface has been observed only in crystals of the 1–2–3 system, where a pronounced guiding effect is expected due to the twin boundaries.

As an illustration we show in Fig. 3 turbulent behaviour observed in an optimally doped YBCO single crystal containing a substantial amount of twinning, see (a). The crystal has a rectangular shape in the ab plane and measures 1 mm along the longest edge. Shown in (b) is the MO-image of flux penetration in an external field of 100 mT applied after zero-field-cooling to 45 K. One sees that the field penetrates predominantly from some large twin boundaries located at the bright core of the lines that make $`45^{}`$-angles with the edges of the crystal. The dark area in the centre of the sample is the region not being penetrated by a 100 mT applied field at this temperature. In (c) the temperature is raised to 67 K, and the MO-image was recorded in the remanent state after first applying 100 mT. The bright ”aura” around the crystal is here the return field of the flux trapped in the central part. Note that this reverse field partly penetrates the sample near the edge. A distinct line can be seen as a dark band going around the crystal just on the inside of the edge. This band is the annihilation zone, which divides the crystal into two opposite magnetic domains. The macroturbulence is here seen as a meandering of the annihilation zone. Adjacent to the zone one can find small areas with increased flux density, see magnified view in (d). Time resolved measurements show that the zone develops in a highly dynamic manner where abrupt redistributions of flux often occur. By adding an external reverse field to the remanent state the annihilation zone is pushed further into the crystal and the dynamical features becomes even more spectacular.

In this crystal the turbulent behavior was observed in the interval 25-75 K. As the temperature increases the dynamics of the flux/antiflux interface becomes increasingly rapid. However, above 75 K there again appears to be no irregular behavior of the interface.

The existence of the instability only in a definite temperature region finds a simple rationalization within our model. At low temperatures, the viscosity increases exponentially, and the characteristic spatial scale $`L`$, in Eq. (9), decreases correspondingly and becomes comparable to or less than the twin-boundary spacing. As a result, the guiding effect is suppressed and the instability disappears. On the other hand, at temperatures close to $`T_c`$ the flux guiding is no longer effective due to thermal activation of the vortices. It is remarkable, and in full support of our model, that in the present heavily twinned crystal the turbulence occurs down to much lower temperatures than found in previous studies of similar crystals with only little twinning .

In conclusion, we have presented a theoretical model which gives a fully consistent description of the macroturbulence phenomenon we observe in HTS materials under certain conditions. This work is supported by INTAS/RFBR and RFBR, grants IR-97–1394 and 00-02-17145, 00-02-18032, and the Research Council of Norway.


format: text

text:
# Dirac particles in a rotating magnetic fieldpublished in J. Phys. A 34 (2001) 1903-1909.

## Acknowledgments

The author is grateful to Professor Guang-jiong Ni for communications and encouragement. This work was supported by the National Natural Science Foundation of China.


text:
# Phase Transition Of The Spin-One Square-Lattice Anisotropic Heisenberg Model

## Abstract

The coupled cluster method (CCM) is applied to the spin-one anisotropic Heisenberg antiferromagnet (HAF) on the square lattice at zero temperature using a new high-order CCM ground-state formalism for general quantum spin number ($`s1/2`$). The results presented constitute the first such application of this new formalism, and they are shown to be among the most accurate results for the ground-state energy and sublattice magnetisation of this model as yet determined. We “track” the solution to the CCM equations at a given level of approximation with respect to the anisotropy parameter, $`\mathrm{\Delta }`$, from the trivial Ising limit ($`\mathrm{\Delta }\mathrm{}`$) down to a critical value $`\mathrm{\Delta }_c`$, at which point this solution terminates. This behaviour is associated with a phase transition in the system, and hence a primary result of these high-order CCM calculations is that they provide an accurate and unbiased (i.e., ab initio) estimation of the position of the quantum phase transition point as a function of the anisotropy parameter. Our result is, namely, that this point occurs at (or slightly below) the isotropic Heisenberg point at $`\mathrm{\Delta }=1`$, for this model.

PACS numbers: 75.10.Jm, 75.50.Ee, 03.65.Ca

The area of quantum magnetic insulating systems at zero temperature has become increasingly well understood for one-dimensional (or quasi-one-dimensional) lattices via the existence of well-known exact solutions such as the Bethe ansatz (BA) ba1 ; ba2  and also via more recent density matrix renormalisation group (DMRG) calculations DMRG1 ; DMRG2 ; DMRG3 . Similarly, zero-temperature quantum Monte Carlo (QMC) calculations qmc3 ; qmc4  have been shown to yield very accurate results for spin-half, two-dimensional systems. QMC techniques, however, are severely limited in their range of application by the presence of the infamous “sign problem,” although we note that for non-frustrated spin systems it is often possible to determine a “sign rule” sign\_rules1 ; sign\_rules2  which completely circumvents this minus-sign problem. By contrast, the technique of quantum many-body theory known as the coupled cluster method (CCM) is neither limited by the presence of frustration nor by the spatial dimensionality of the lattice ccm2 ; ccm7 ; ccm4 ; ccm5 ; ccmextra ; ccm6 . Indeed, a great strength of this method is that it is able to determine with great accuracy the positions of quantum phase transition points of quantum systems within an ab initio, and thus essentially unbiased, framework. These quantum phase transitions typically arise as some parameter within the Hamiltonian is varied, thus driving the system from one phase into another.

The behaviour of the quantum Heisenberg antiferromagnet (HAF) on the linear chain is highly dependent on whether the quantum spin number, $`s`$, has an integer or a half-odd-integer value. Indeed, Haldane haldane  first predicted that the spin-one isotropic HAF would contain an excitation energy gap, and this prediction has subsequently and conclusively been shown numerically to be correct by exact diagonalisations of short chains parkinson  and by more recent and extremely accurate DMRG calculations DMRG3 . We note that this is in stark contrast to the exact BA solution ba1 ; ba2  of the spin-half isotropic HAF for the linear chain, which contains no such gap. However, such exact diagonalisations and DMRG calculations are considerably more difficult for systems of higher spatial dimensionality. In this article, we focus on the specific case of the spin-one anisotropic HAF on the square lattice at zero temperature using the CCM, especially in relation to its quantum phase transition, by applying a high-order CCM formalism for general quantum spin number for the first time.

The spin-one anisotropic HAF is given by,

$$H=\underset{i,j}{}\left\{s_i^xs_j^x+s_i^ys_j^y+\mathrm{\Delta }s_i^zs_j^z\right\},$$
(1)
where the symbol $`i,j`$ indicates nearest-neighbour bonds on the square lattice, and where each bond is counted once and once only. For $`\mathrm{\Delta }\stackrel{>}{_{}}1`$, the ground state of this model contains non-zero Néel-like order. The precise value of the phase transition point, at which this ordering breaks down is not exactly known although it is believed to be at (or near to) the Heisenberg point at $`\mathrm{\Delta }=1`$. Very accurate results for the values of various ground-state properties of this model have been obtained via spin-wave theory (SWT) Hamer  and cumulant series expansions Zheng . Both sets of results indicate that approximately $`80\%`$ of the classical ordering remains for the spin-one, square-lattice HAF model at $`\mathrm{\Delta }=1`$. We note however that cumulant series expansions make an explicit assumption that the position of the phase transition point is exactly at $`\mathrm{\Delta }=1`$ in order to perform Padé resummations of the otherwise divergent perturbation series. Also, conventional SWT essentially appears to “build in” a phase transition point at $`\mathrm{\Delta }=1`$ for any lattice, at which point the excitation spectrum becomes soft. This is, of course, now known to be incorrect DMRG3 ; haldane ; parkinson  for the spin-one linear-chain HAF model, as noted above.

We now briefly describe the general CCM formalism, and the interested reader is referred to Refs. ccm2 ; ccm7 ; ccm4 ; ccm5 ; ccmextra ; ccm6  for further details. The exact ket and bra ground-state energy eigenvectors, $`|\mathrm{\Psi }`$ and $`\stackrel{~}{\mathrm{\Psi }}|`$, of a many-body system described by a Hamiltonian $`H`$ are parametrised within the single-reference CCM as follows:

$`|\mathrm{\Psi }=\mathrm{e}^S|\mathrm{\Phi }`$ ; $`S={\displaystyle \underset{I0}{}}𝒮_IC_I^+,`$
$`\stackrel{~}{\mathrm{\Psi }}|=\mathrm{\Phi }|\stackrel{~}{S}\mathrm{e}^S`$ ; $`\stackrel{~}{S}=1+{\displaystyle \underset{I0}{}}\stackrel{~}{𝒮}_IC_I^{}.`$ (2)
The single model or reference state $`|\mathrm{\Phi }`$ is required to have the property of being a cyclic vector with respect to two well-defined Abelian subalgebras of multi-configurational creation operators $`\{C_I^+\}`$ and their Hermitian-adjoint destruction counterparts $`\{C_I^{}(C_I^+)^{}\}`$. Thus, $`|\mathrm{\Phi }`$ plays the role of a vacuum state with respect to a suitable set of (mutually commuting) many-body creation operators $`\{C_I^+\}`$. Note that $`C_I^{}|\mathrm{\Phi }=0,I0`$, and that $`C_0^{}1`$, the identity operator.

We note that the exponentiated form of the ground-state CCM parametrisation of Eq. (2) ensures the correct counting of the independent and excited correlated many-body clusters with respect to $`|\mathrm{\Phi }`$ which are present in the exact ground state $`|\mathrm{\Psi }`$. It also ensures the exact incorporation of the Goldstone linked-cluster theorem, which itself guarantees the size-extensivity of all relevant extensive physical quantities.

The CCM equations are defined by the following coupled set of equations,

$`\mathrm{\Phi }|C_I^{}\mathrm{e}^SH\mathrm{e}^S|\mathrm{\Phi }`$ $`=`$ $`0,I0;`$ (3)
$`\mathrm{\Phi }|\stackrel{~}{S}\mathrm{e}^S[H,C_I^+]\mathrm{e}^S|\mathrm{\Phi }`$ $`=`$ $`0,I0.`$ (4)
We furthermore note that the ground-state energy at the stationary point has the simple form

$$E_g=E_g(\{𝒮_I\})=\mathrm{\Phi }|\mathrm{e}^SH\mathrm{e}^S|\mathrm{\Phi }.$$
(5)
It is important to realize that this (bi-)variational formulation does not lead to an upper bound for $`E_g`$ when the summations for $`S`$ and $`\stackrel{~}{S}`$ in Eq. (2) are truncated, due to the lack of exact Hermiticity when such approximations are made. However, it is clear that the important Hellmann-Feynman theorem is still preserved in all such approximations.

The CCM formalism is exact in the limit of inclusion of all possible multi-spin cluster correlations for $`S`$ and $`\stackrel{~}{S}`$, although in any real application this is usually impossible to achieve. It is therefore necessary to utilise various approximation schemes within $`S`$ and $`\stackrel{~}{S}`$. The three most commonly employed schemes previously utilised have been: (1) the SUB$`n`$ scheme, in which all correlations involving only $`n`$ or fewer spins are retained, but no further restriction is made concerning their spatial separation on the lattice; (2) the SUB$`n`$-$`m`$ sub-approximation, in which all SUB$`n`$ correlations spanning a range of no more than $`m`$ adjacent lattice sites are retained; and (3) the localised LSUB$`m`$ scheme, in which all multi-spin correlations over distinct locales on the lattice defined by $`m`$ or fewer contiguous sites are retained. We note that in order to carry out such LSUB$`m`$ or SUB$`n`$-$`m`$ calculations to high order in the truncation indexes $`n`$ and $`m`$ we must rely on computer-algebra techniques in order to generate the corresponding sets of coupled equations. These computational techniques are based on the fact that the similarity transformed Hamiltonian in Eqs. (3-5) may be written purely in terms of creation operators acting on the model state, so that evaluation of Eq. (3) becomes purely a matter of pattern matching. The interested reader is referred to Refs. ccm4  and ccm6  for a full explanation of how this is achieved for spin-half models. We note that the generalisation of this procedure to the case of arbitrary quantum spin number, $`s`$, is relatively straightforward, although detailed, and will be described fully elsewhere ccm8 .

We now apply the CCM formalism outlined above to the specific case of the anisotropic HAF, and we choose the Néel state, in which the spins lie along the $`z`$-axis, to be the model state. Furthermore, we perform a rotation of the local axes of the up-pointing spins by 180 about the $`y`$ axis, so that spins on both sublattices may be treated equivalently. The (canonical) transformation is described by,

$$s^xs^x,s^ys^y,s^zs^z.$$
(6)
The model state now appears $`mathematically`$ to consist of purely down-pointing spins in these rotated local axes. In terms of the spin raising and lowering operators $`s_k^\pm s_k^x\pm \mathrm{i}s_k^y`$ the Hamiltonian may be written in these local axes as,

$$H=\frac{1}{4}\underset{i,j}{}\left\{s_i^+s_j^++s_i^{}s_j^{}+2\mathrm{\Delta }s_i^zs_j^z\right\},$$
(7)
where $`i`$ and $`j`$ runs over all nearest neighbours, although each nearest-neighbour bond is counted once and once only. Furthermore, the sublattice magnetisation, $`M`$, (after rotation of the local spin axes) is given by,

$$M=\frac{1}{N}\underset{i=1}{\overset{N}{}}s_i^z.$$
(8)
In the limit $`\mathrm{\Delta }\mathrm{}`$ we note that the model state is the exact ground eigenstate of the Hamiltonian of Eq. (7). Hence, all of the CCM correlation coefficients are zero at this point. We may track this solution for decreasing values of $`\mathrm{\Delta }`$ until we reach a critical point, $`\mathrm{\Delta }_c`$, at which point the real solution to our CCM equations for the LSUB$`m`$ and SUB$`m`$-$`m`$ approximation schemes terminates. This is associated with a phase transition in the real system ccm4 ; ccm5 ; ccmextra ; ccm6 . Similar behaviour was also seen for this model for the CCM SUB2 approximation for this model (see Ref. ccm7 ), and once more this was associated with a phase transition of the real system.

We note that the “raw” results for the ground-state energy, sublattice magnetisation, and critical points, $`\mathrm{\Delta }_c`$, have been obtained using both the LSUB$`m`$ and SUB$`m`$-$`m`$ approximation schemes. For both approximation schemes we may extrapolate these results as a function of $`m`$ in the limit $`m\mathrm{}`$ in order to obtain even better results. We note that SUB2-$`m`$ calculations and the full SUB2 calculation have also previously been performed ccm7  for the anisotropic square-lattice HAF with general quantum spin number. It was therefore possible to determine accurately the manner in which SUB2-$`m`$ results scale as a function of $`m`$, namely, that as $`m\mathrm{}`$ the ground-state energy and the critical points, $`E_g(m)/N`$ and $`\mathrm{\Delta }_c(m)`$ respectively, scale linearly with $`1/m^2`$, and the sublattice magnetisation, $`M(m)`$, scales linearly with $`1/m`$. By analogy, we utilise a similar procedure here in order to extrapolate raw LSUB$`m`$ and SUB$`m`$-$`m`$ results. However, as we may utilise LSUB$`m`$ and SUB$`m`$-$`m`$ results for $`m=\{2,4,6\}`$ only, a quadratic function is fitted to these data in order to obtain the best possible results. A full and comprehensive explanation of the extrapolation process of CCM LSUB$`m`$ expectation values is given in Ref. ccm6  for the spin-half XXZ model for a variety of lattices, and the interested reader is referred to this article for further information.

CCM results for the position of the critical point using the LSUB$`m`$ and SUB$`m`$-$`m`$ approximation schemes are presented in Tables 1 and 2. The extrapolated result of $`\mathrm{\Delta }_c=0.98`$ for both the LSUB$`m`$ and for the SUB$`m`$-$`m`$ approximation scheme indicates that the phase transition point is at (or perhaps slightly below) the isotropic Heisenberg point ($`\mathrm{\Delta }=1`$). The strength of the CCM is that it can provide such an accurate value for the position of the quantum phase transition point as a function of the anisotropy within an ab initio, and thus fully unbiased, framework.

Figure 1 and Tables 1 and 2 indicate that the results for the ground-state energy of the spin-one HAF converge extremely quickly with $`m`$ for both the LSUB$`m`$ and SUB$`m`$-$`m`$ schemes. Indeed, even the “raw” unextrapolated results provide excellent estimates of the ground-state energy, although the extrapolated results of $`E_g/N=2.3292`$ for the LSUB$`m`$ approximation scheme and $`E_g/N=2.3291`$ for SUB$`m`$-$`m`$ approximate scheme are certainly even more accurate. By way of comparison, we note that third-order SWT Hamer  gives a result of $`E_g/N=2.3282`$ and cumulant series expansions Zheng  a result of $`E_g/N=2.3279(2)`$.

CCM results for the sublattice magnetisation are found to be similarly well converged as a function of the truncation index $`m`$ for both the LSUB$`m`$ and SUB$`m`$-$`m`$ schemes, as indicated in Fig. 2. The extrapolated results of $`M=0.8049`$ and $`M=0.8067`$ for the LSUB$`m`$ and SUB$`m`$-$`m`$ schemes, respectively, are again in excellent agreement with the results of third-order SWT Hamer  and cumulant series expansions Zheng , which give values of $`M=0.8043`$ and $`M=0.8039(4)`$ respectively.

In this article it has been shown that the coupled cluster method (CCM) provides quantitatively accurate results for the ground-state properties of the spin-one square-lattice isotropic HAF by comparison with results of third-order spin-wave theory (SWT) and cumulant series expansions. Furthermore, we note that the results presented in this article constitute the first application of the CCM using a new high-order CCM ground-state formalism for general ($`s1/2`$) quantum spin number. (A fuller account of this new high-order formalism will be published elsewhere ccm8 .) The best estimates of the ground-state energy of this model were found to be $`E_g/N=2.3292`$ using the LSUB$`m`$ approximation scheme and $`E_g/N=2.3291`$ using the SUB$`m`$-$`m`$ approximation scheme via heuristic extrapolation to the (exact) limit $`m\mathrm{}`$. The best estimates of the sublattice magnetisation are $`M=0.8049`$ using the LSUB$`m`$ approximation scheme and $`M=0.8067`$ using the SUB$`m`$-$`m`$ approximation scheme, again via extrapolation to the limit $`m\mathrm{}`$. The most important result of the CCM calculations presented here for this model is the prediction that the phase transition point of the spin-one square-lattice anisotropic HAF is at (or slightly below) the isotropic Heisenberg point. A strength of the method is that this prediction may be made using an unbiased (i.e., ab initio) treatment.


text:
# Untitled Document

TEXT WITHDRAWN


text:
# Incorporating next-to-leading order matrix elements for hadronic diboson production in showering event generators

## I INTRODUCTION

As high energy collider physics moves towards the TeV frontier, the simulation of collision events at higher orders becomes increasingly important. For the case of hadronic diboson production, higher order corrections are largest in the regions where the physics is most interesting. Event rate enhancements at next-to-leading order in QCD (NLO) have a spoiling effect on the confidence limits for anomalous triple gauge boson couplings, wash out the radiation zero in $`W\gamma `$ production  and approximate radiation zero in $`WZ`$ production , and contribute to the background rate for new physics channels including heavy Higgs boson searches and Supersymmetry signatures.

Until recently, general purpose event generators have included leading order (LO) matrix elements only. An exclusive (meaning all of the parton emissions are explicitly enumerated in the event record) all orders description of multiple emissions is provided by the parton shower approach  which is valid in the regions of soft or collinear emissions, but is not accurate for well-separated partons. The shower occurs with unit probability and so does not alter the cross section which remains leading order. Two of the general purpose showering Monte Carlo programs, HERWIG and PYTHIA, have implemented matrix element corrections  to the parton shower algorithms using the first order tree level matrix element (the one-loop virtual graphs are not included). The corrections are successful in improving the agreement with data in distribution shapes, especially at high transverse momentum. The normalization, however, is everywhere leading order and the overall event rate is usually obtained by normalizing to the data or applying a constant factor obtained from NLO calculations.

In contrast, next-to-leading (fixed) order Monte Carlo integration programs such as the NLO calculations available for $`WZ`$ production  provide a first order prediction of the cross section and a good description of the hardest emission. However, these programs frequently generate events with negative weights in certain phase space regions, do not address multiple emission, and give a poor (indeed sometimes unphysical) description in the region of very small jet transverse momentum where perturbative QCD is unreliable. In this small transverse momentum region, the effects of multiple emissions can be resummed (first proposed by ) giving a better prediction of inclusive An example of an inclusive quantity is the color singlet system transverse recoil against all parton emissions ($`P_{\mathrm{W}Z}^T`$ for the case of $`WZ`$ production). Inclusive means the contribution from individual emissions is not available, only the contribution summed over all emissions is known. A LO cross section is an inclusive prediction of the event rate. After applying a parton shower algorithm to a LO event, the description becomes exclusive because the enumeration of each parton emission is contained in the event record. quantities. However, the standard resummation formalism does not describe the exclusive set of partons in the event. In Ref.  a resummed treatment of hadronic $`Z`$-pair production $`ppZZ`$, is given and a comparison is made with the PYTHIA showering event generator. In the region of low $`Z`$-pair transverse momentum $`P_{\mathrm{Z}Z}^T`$, the resummed spectrum agrees reasonably well with the results from the parton shower approach.

The combination of NLO matrix elements with parton showering algorithms is a natural way to proceed and would provide an exclusive Monte Carlo description of collision events. This goal has two primary challenges: (1) A probabilistic interpretation of NLO matrix elements is difficult to accomplish on account of the delicate cancellations which occur between the virtual and real emission graphs. This generally results in event weights which may be both positive or negative. (2) The parton shower algorithm can generate emissions which are harder than the first (fixed order) emission resulting in a double counting of hard emissions. Both of these drawbacks are accentuated in the region of small inclusive jet transverse momentum.

Ref.  addresses the first challenge for the special case of hadronic $`WZ`$ production with leptonic decays. For each event, a volume of low transverse momentum ($`P_{}^T`$) parton emissions (defined by a $`P_{}^T`$ cutoff and assumed to be unobservable by a detector) is integrated numerically keeping the observable quantities (the charged lepton momenta and missing transverse momentum) fixed. The integration volume is generated by allowing the missing transverse momentum to arise from any combination of invisible neutrino and low $`P_{}^T`$ emission. Such an integration ensures one doesn’t “look into” the region of small inclusive jet transverse recoil $`P_{\mathrm{W}Z}^T`$, where the fixed order calculation is unreliable. The resultant events can be interpreted probabilistically (event weights are positive definite) and the NLO predictions for observable distributions are preserved. The drawback is that by integrating the volume of jets the specific jet configuration is lost—determining which configuration is most probable brings back the problem of negative weights. This, and the discontinuity caused by the $`P_{}^T`$ cutoff, makes the subsequent application of a parton shower algorithm difficult.

In the present study the positive features of the Ref.  algorithm are combined with the matrix element corrections approach. The phase space is divided into two complementary regions: the region accessible to the parton shower (“parton shower region”) and the hard central region which is not accessible to the parton shower (“Dead Zone region”). The treatment of the Dead Zone region in this study follows Ref. , wherein matrix element corrections to single vector boson production are applied within the framework of the HERWIG program. The extension to diboson production is trivial, the vector boson mass and momentum in Ref.  is everywhere replaced by the diboson system mass and momentum. The parton shower region serves as a phase space volume over which emissions can be integrated in a treatment that closely follows Ref. . This provides a means of including information from the full NLO matrix element (including one-loop graphs) such that the normalization is NLO and the distributions reproduce the NLO shapes in the regions where the NLO prediction is expected to be reliable. In the low inclusive jet transverse momentum region, where multiple emissions are important and the fixed order prediction becomes unreliable, the distribution shapes are dominated by the all orders parton shower provided by HERWIG (though the normalization remains NLO). In this way events generated with the algorithm described here provide an exclusive description of the jet structure while giving a physical result across the full jet transverse momentum spectrum and providing NLO normalization.

The method described in this paper applies to any hadronic color singlet production process. The specific process of $`ppW^+Z`$ production with vector boson decays to electron/muon-type leptons at the Large Hadron Collider (LHC) energy of 14 TeV is chosen to illustrate the method. The input parameters for this study are $`\alpha _{EM}(M_Z)=\frac{1}{128}`$, $`\mathrm{sin}^2\theta _W=0.23`$, $`\alpha _s(M_Z)=0.1116`$, $`M_W=80.396`$ GeV, $`M_Z=91.187`$ GeV, Cabibbo angle $`\mathrm{cos}\theta _{\mathrm{C}abibbo}=0.975`$ with no 3rd generation mixing, and the factorization scale $`Q`$ is taken to be the hard process mass $`\sqrt{\widehat{s}}`$. Vector boson branching ratios to leptons are taken as Br($`Zl^+l^{}`$) = 3.36%, Br($`W^\pm l^\pm \nu `$)=10.8%. The CTEQ4M  structure functions are used with the $`b`$ quark contribution taken as zero. Matrix elements (LO and NLO) are evaluated with the Baur-Han-Ohnemus (BHO) package  and HERWIG is employed for the parton cascade, hadronization, and decays.

In the section which follows, the implementation of matrix element corrections for diboson production in the HERWIG style are presented and a comparison is made between LO with parton shower, NLO, and matrix element corrected distributions. The method for incorporating the full NLO matrix element in the HERWIG showering event generator is presented in Sec. III and is compared with matrix element corrected and NLO distributions. Conclusions are given in the final section.

## II MATRIX ELEMENT CORRECTIONS

The parton level Born process for hadronic $`WZ`$ production is $`q\overline{q}^{}WZ`$. The kinematics describing this process are given the label ‘n-body’ and the label ‘(n+1)-body’ is reserved for the first order tree processes $`q\overline{q}^{}WZg`$, $`qgWZq^{}`$, and $`g\overline{q}WZ\overline{q}^{}`$.

Given the kinematics of an n-body event as input, a showering Monte Carlo program uses backwards evolution to trace the hard scattering partons back from the hard vertex into the original incoming hadrons while keeping track of the distribution of emitted partons (refer to Ref.  for a good review). For the case of the HERWIG program, the phase space accessible to the parton shower (the parton shower region) is a subset of the full (n+1)-body phase space. The boundary of the Dead Zone region inaccessible to the HERWIG parton shower is defined by

$$\frac{7\sqrt{17}}{2}M^2<\widehat{s}<s$$
(1)
and

$$|\mathrm{cos}\widehat{\theta }_j|<1\frac{3\sqrt{1+8\frac{M^2}{\widehat{s}}}}{\frac{\widehat{s}}{M^2}1}$$
(2)
where $`M`$ is the mass of the color singlet object ($`WZ`$ system in this case), $`\sqrt{\widehat{s}}`$ is the mass of the color singlet and parton emission system ($`WZj`$), $`\sqrt{s}`$ is the machine energy, and $`\widehat{\theta }_j`$ is the angle of the emission with respect to the beam in the center of mass frame. The partition of phase space in the $`\widehat{s}/M^2`$ vs. $`\mathrm{cos}\widehat{\theta }_j`$ plane is shown in Figure 2. The parton shower and Dead Zone regions are precisely complementary, so there can be no double counting when sampling the regions separately. The Dead Zone is the region of hard central emissions and can be populated by sampling the first order tree level matrix element within the boundaries defined in Eqs. 1 and 2. Adding the Dead Zone phase space in this manner to the showering event generator simulation is termed ‘hard matrix element corrections’. For single boson production the improvement in the boson recoil transverse momentum spectrum  is remarkable, particularly at high $`P_{}^T`$.

‘Soft matrix element corrections’  refers to correcting emissions generated by the parton shower algorithm using the first order tree level matrix element. Since the HERWIG parton shower ordering variable is an energy-weighted angle which does not necessarily imply ordering in transverse momentum, it is necessary to apply the correction to every emission which has the highest $`P_{}^T`$ so far. The soft corrections are described for the single boson production case in Ref. and were noted to have a small effect in comparison to the hard matrix element corrections for single boson production. The main advantage of the soft corrections is in providing for a smooth transition between the parton shower and Dead Zone regions—though in practice it is difficult to find observables sensitive to this transition. The primary difference with regards to soft corrections for diboson production as compared to the single boson case are the real emission graphs involving a gluon anchored to an internal line, such as the one shown in Figure 1. These graphs are absent for the single boson case and are in no way approximated by the parton shower which traces only external lines. In contrast to the other emission graphs, these graph do not have any collinear singularities. This means that their relative contribution to the cross section is largest in the central regions (since other emission graphs dominate in the collinear regions). Thus, the relative effect of these graphs in the parton shower region is small, and it is more important to account for them in the Dead Zone region where the physics of these graphs may become important.

The present study uses the soft matrix element corrections already implemented in HERWIG for single boson production. Since the emission cross section enters only as a ratio to the Born cross section, the HERWIG implementation accounts for the first order emission from an external line of Feynman graphs for color singlet production to good approximation, but does not account for emissions from an internal line. Implementing the soft corrections for the specific case of $`WZ`$ production would require intrusive changes to the HERWIG parton shower code. By accepting this approximation, the HERWIG code is kept unmodified and modular. Nevertheless it is recognized that explicit diboson soft matrix element corrections would be a useful extension to this study.

The HERWIG soft matrix element corrections (which are approximate for the diboson case) have very little impact on the distributions presented in this paper. It is expected that the exact soft matrix element corrections (i.e. explicitly including all diboson graphs) would not change the situation in an appreciable way.

### A  Matrix element corrected results

LO and matrix element corrected distributions for hadronic $`WZ`$ production are obtained by implementing the Born level n-body process and tree level (n+1)-body process (restricted to the Dead Zone) as external processes in HERWIG using the BHO matrix elements. For diboson production, the distribution corresponding to the single vector boson production transverse recoil studied in Ref.  is the diboson system recoil $`P_{\mathrm{W}Z}^T`$. A comparison to the NLO result for this distribution is shown in Figure 3 for the LO with parton shower, first order tree level in the Dead Zone, and matrix element corrected results. Away from the low $`P_{}^T`$ region, the shape of the matrix element corrected result agrees well with the NLO prediction from the stand-alone BHO integration package, demonstrating the effectiveness of hard matrix element corrections. There is no contribution from the Dead Zone in the very low $`P_{\mathrm{W}Z}^T`$ region, so the Born and matrix element corrected results coincide in this region (Figure 3, inset). The soft matrix element corrected LO result is indistinguishable from the LO curve and is not shown.

It is more common in the literature to study the $`P_\mathrm{Z}^T`$ distribution, for example in probing the triple gauge boson vertex for anomalous couplings. This distribution is also sensitive to hard emissions in the Dead Zone, as shown in Figure 4 where, as for the $`P_{\mathrm{W}V}^T`$ spectrum, the shape of the matrix element corrected $`P_\mathrm{Z}^T`$ distribution is a good approximation to the NLO result.

Drawbacks in the matrix element corrections approach and fixed order calculation are evident in Figure 3. The NLO calculation is not reliable in the very low $`P_{}^T`$ region shown in the inset of Figure 3. The NLO prediction shows a dip just below 5 GeV and a very large contribution at zero. These features are not physical and reflect the specific choices of the regularization scheme used in the matrix elements. A different choice of scheme or scheme parameters would yield a very different result in this region. The BHO matrix elements employ the two parameter phase space slicing method  as regularization scheme. The specific choice of cutoff parameters used for Figure 3 are large ($`\delta _\mathrm{s}=0.05,\delta _\mathrm{c}=0.01`$ see Ref. for a definition of the cutoffs). A smaller choice of parameters would move events out of the zero histogram bin of the Figure 3 inset towards regions of higher $`P_{}^T`$, until eventually the zero bin would become negative. The use of the subtraction method  for the $`P_{\mathrm{W}Z}^T`$ distribution in the inset of Figure 3 also yields an unphysical negative first histogram bin. This effect serves as a reminder that fixed order perturbation theory is not well suited for predictions of small inclusive $`P_{}^T`$.

A disadvantage of the matrix element corrected distributions is the LO normalization in the parton shower region. Since the total matrix element corrected cross section is the sum of the Born level and the integrated Dead Zone cross sections, it depends explicitly on the Dead Zone boundary defined by Eqs. 1 and 2. A change in the boundary would result in a change in the total cross section. The use of NLO calculations everywhere would remove this dependence and provide a more accurate prediction of the total cross section, particularly for multi-TeV colliders where the contribution from the gluon (anti-)quark initial state is substantial. A second drawback of the matrix element corrections is the inability to model NLO shape effects in the parton shower region. The approximate radiation zero, which is partially filled in by NLO corrections, is an example of a distribution which is sensitive to these effects. Observation of this as yet unseen phenomenon will be a goal of future hadron colliders. This distribution is shown in Figure 5, together with the ratio (or $`k`$-factors) of the matrix element corrected and Born distributions with respect to the NLO distributions. The $`k`$-factors are largest in the central radiation zero region and the matrix element corrections are successful in greatly reducing this effect. However, a small variation in $`k`$-factors in the zero rapidity separation region remains after applying the matrix element corrections.

## III INCORPORATING THE FULL NLO MATRIX ELEMENT

In the Dead Zone the square of the real emission graphs ($`_{\mathrm{r}ealemission}^2`$) is the only contribution to the cross section and the divergent regions of the real emission phase space are excluded. Outside of the Dead Zone the parton shower is better suited at predicting event shapes and exclusive distributions, though it provides no information on the normalization which in showering event generators is normally derived from the LO cross section only. In this region it is preferable to obtain the normalization from the full NLO cross section while preserving the exclusive event structure from the parton shower.

In the parton shower region the square of the Born graphs, the interference of the Born graphs with the one-loop graphs, and the square of the real emissions graphs contribute at NLO,

$$_{\mathrm{N}LO}^2=_{\mathrm{B}orn}^2+_{\mathrm{B}orn}_{\mathrm{o}neloop}+_{\mathrm{r}ealemission}^2.$$
(3)
The first two terms on the right hand side of Eq. 3 are described by n-body kinematics and occupy a point at $`(\widehat{s}/M^2,\mathrm{cos}\widehat{\theta }_j)=(1,0)`$ of Figure 2. The last term is described by (n+1)-body kinematics and may populate the entire $`\widehat{s}/M^2`$ vs. $`\mathrm{cos}\widehat{\theta }_j`$ plane. A regularization scheme is necessary to handle the soft, collinear, and ultraviolet divergences which appear when any of the first order contributions are treated alone. The BHO matrix elements use the phase space slicing method  (PSS) wherein a small portion (denoted on Figure 2) of the (n+1)-body phase space (defined by soft and collinear cutoff parameters $`\delta _\mathrm{s},\delta _\mathrm{c}`$) is included in the n-body contribution by means of the soft gluon and leading pole approximations.

One way of achieving NLO normalization while gaining sensitivity to the effect of NLO corrections on the $`WZ`$ system configuration is to assign an event weight for the hard process which is the average NLO cross section (including virtual graphs) across the entire region accessible to the parton shower for a particular choice of the $`WZ`$ system configuration. Since all of the divergent regions are contained within the parton shower region, this number is well behaved and positive definite—but is difficult to calculate analytically.

This average cross section can be obtained by means of a 2nd stage integration as outlined in Ref. . It is also noted that it is not necessary to obtain this number with a high degree of accuracy—so long as the average of many such integrations (i.e. events) converges to the correct cross section, which is automatic using the Monte Carlo method.

For the present algorithm the jet volume has a simple definition (Eqs. 1 and 2) with the convenient feature that the phase space slicing method jet volume is a subset of this volume. This means the phase space slicing method can be used to achieve a large portion of the jet volume integration and that the result is not sensitive to the specific choice of PSS cutoff parameters which specify the PSS boundary. As such, the maximal PSS cutoffs (which remain within the range of validity of the soft gluon and leading pole approximations) may be used ($`\delta _\mathrm{s}=0.05,\delta _\mathrm{c}=0.01`$ have been used, see Ref. for a definition of the cutoffs). This is a significant advantage because the subset occupied by the PSS volume is the region where delicate cancellations between the n-body and (n+1)-body graphs occur and also where the cross section is peaked.

The 2nd stage integration, which is performed once for each event in the parton shower region, is accomplished as follows. A $`WZ`$ event with n-body kinematics is generated (8 degrees of freedom specify the n-body configuration assuming the vector bosons are taken on shell: e.g. the $`WZ`$ system mass $`M_{\mathrm{W}Z}`$, $`WZ`$ system rapidity $`y_{\mathrm{W}Z}`$, the 2 angles defining the vector boson production, and 2 angles for each of the vector boson decays) and defines the $`WZ`$ system configuration and boost. For a LO showering Monte Carlo event, the Born level cross section would be evaluated and constitute the event weight. Instead, the aim is to evaluate the average cross section of the phase space into which the parton shower may evolve the event.

The NLO n-body cross section is calculated using the BHO matrix element which requires a further 4 degrees of freedom to specify the soft and collinear corrections arising in the PSS method. Thus by integrating over these 4 degrees of freedom the average n-body weight $`d\sigma ^{\mathrm{n}\mathrm{body}}`$ is obtained and accounts for the region inside the PSS boundary (on the left side of the dashed line) in Figure 2. This number may be negative. For the integration of the remainder of the parton shower region an emission is explicitly sampled using 3 additional degrees of freedom (e.g. the $`WZj`$ system mass and the 2 angles specifying the emission direction). These 11 degrees of freedom (3 plus 8 from the n-body configuration) specify the (n+1)-body kinematics. It is necessary to use the appropriate overall boost, so as to span the identical phase space the HERWIG parton shower is able to reach from the n-body kinematics. HERWIG applies a conformal boost to the hard process which keeps the hard process mass and rapidity fixed. For the (n+1)-body case the $`WZj`$ system rapidity is thus

$$y_{(\mathrm{n}+1)\mathrm{body}}=\mathrm{ln}\left(\frac{\sqrt{M_{\mathrm{W}V}^{}{}_{}{}^{2}+P_{\mathrm{W}Z}^{T}{}_{}{}^{2}}}{E_{\mathrm{W}Z}^{}{}_{}{}^{}+P_{\mathrm{W}Z}^{z}{}_{}{}^{}}\right)+y_{\mathrm{n}\mathrm{body}}$$
(4)
where $`y_{\mathrm{n}\mathrm{body}}`$ is identified as the $`WZ`$ system rapidity and $`P_{\mathrm{W}Z}^{}{}_{}{}^{}`$ is the $`WZ`$ system 4-vector in the $`WZj`$ frame. The Jacobian for the transformation $`y_{\mathrm{n}\mathrm{body}}y_{(\mathrm{n}+1)\mathrm{body}}`$ is unity. By sampling the 3 degrees of freedom which specify the expansion from n-body to (n+1)-body kinematics and evaluating the NLO matrix element each time (always a positive number), an integration over the remainder of the parton shower region is achieved giving the average (n+1)-body weight $`d\sigma ^{(\mathrm{n}+1)\mathrm{body}}`$.

The cross section for the specific choice of $`WZ`$ configuration averaged over the parton shower region at NLO is

$$d\sigma ^{\mathrm{p}artonshowerregion}=d\sigma ^{\mathrm{n}\mathrm{body}}+d\sigma ^{(\mathrm{n}+1)\mathrm{body}}$$
(5)
and constitutes the event weight. The input to the parton shower algorithm is the n-body kinematic configuration. The parton shower algorithm (with guidance from the soft matrix element corrections) determines which exclusive emission structure is chosen for the event.

Event generation for the NLO $`WZ`$ process consists of generating a sample of events in the parton shower region using the 2nd stage integration method discussed above and a second sample of events in the Dead Zone, as already discussed in Sec. II. The event weights for both types of events are positive numbers (unlike for the NLO integration packages which produces events with both positive and negative weights) so the usual event generator strategy of hit-and-miss may be used to obtain unweighted event samples.

### A Results

The cross sections for the regions of interest are enumerated in Table I for $`W^+Z`$ production at LHC with decays to electron/muon-type leptons. The total Born level cross section is almost a ($`k`$-)factor 1.5 down from the NLO result. The Dead Zone accounts for about 8% of the total cross section. The sum of the integrated parton shower region cross section and the Dead Zone cross section agrees with the total NLO cross section which serves as a cross check that the parton shower region contribution is correctly evaluated with the 2nd stage integration.

In Figure 6 the transverse recoil of the $`WZ`$ system is shown for this 2nd stage integration method (which also includes matrix element corrections). The normalization is the same as the NLO normalization by construction. The NLO and matrix element corrected results from Figure 3 are superimposed. In the region of moderate to high $`P_{\mathrm{W}Z}^T`$ the result agrees well with the NLO prediction. In the region of smaller $`P_{\mathrm{W}Z}^T`$, events with low transverse momentum receive extra recoil from multiple emissions and are pushed towards slightly higher $`P_{\mathrm{W}Z}^T`$. At very small $`P_{\mathrm{W}Z}^T`$ (inset) the shape agrees very well with the parton shower distribution and goes smoothly to zero with no dependence on the NLO regularization scheme. The ratio to the fixed NLO distribution is slightly flatter (Fig. 6, bottom) than for the matrix element corrected distribution.

A second distribution which is sensitive to the region of small inclusive jet $`P_{}^T`$ is the azimuthal separation of the vector bosons shown in Figure 7. As the separation approaches $`\pi `$, the inclusive jet $`P_{}^T`$ goes to zero and the fixed order prediction becomes unreliable, as evidenced by the sharp increase in differential cross section very near $`\varphi (W)\varphi (Z)=\pi `$. The 2nd stage integration and matrix element corrected results do not suffer from this discontinuity and give a smooth, physical prediction across the full range which is higher than the NLO result away from $`\varphi (W)\varphi (Z)=\pi `$. This is an example of a distribution where the NLO prediction is not accurate.

The $`Z`$ transverse momentum is shown in Figure 8 for the 2nd stage integration method, the matrix element corrections, and the NLO calculation. The 2nd stage integration is in good agreement with the NLO prediction throughout, while the matrix element corrections suffer from normalization problems. The NLO calculation is accurate in the low $`P_\mathrm{Z}^T`$ region for this observable since an integration over a large range of inclusive jet $`P_{}^T`$ is automatic for the low $`P_\mathrm{Z}^T`$ bins.

The rapidity separation of the vector bosons is shown in Figures 9 for the 2nd stage integration method, the matrix element corrections, and the NLO calculation. The 2nd stage integration is in good agreement with the NLO prediction throughout, with a very slight enhancement visible in the radiation zero region due to multiple emissions. The ratio of the 2nd stage integration distribution to the NLO result is unity and almost constant throughout, whereas the matrix element corrected ratio is about 1.3 with a small variation in the zero separation region. The 2nd stage integration has accounted for the effects of the first order corrections on the $`WZ`$ system configuration (including the separation in rapidity) and so this small variation is not present in that distribution.

### B  Complications due to negative weights

For all of the 2nd stage integration method distributions (Figures 6-9) the integration over the parton shower region (performed once per event) is accomplished with 50 samples of the 3 degrees of freedom which define the expansion to (n+1)-body phase space and 10 samples of the 4 degrees of freedom which define the n-body PSS corrections. The 2nd stage integration increases the computer time required to generate one weighted event with the HERWIG showering event generator by a factor 1.3.<sup>§</sup><sup>§</sup>§ The computation time will increase significantly when generating unweighted events.  This is due mostly to the integration, but also to the added number of mathematical operations necessary to evaluate the NLO matrix element as compared to the LO one.

Since the n-body contribution can be negative, the total event weight (Eq. 5) can be within statistical precision of zero if the integration statistics (number of samples) are small enough and the cancellations (between n-body and (n+1)-body contributions) are large enough. This means that negative events can very occasionally occur (a discussion of 2nd stage integration negative weighted events is given in Ref.  The lepton correlation effect discussed in Ref.  also occurs in the present study, but the effect is small. ). Negative events of this type are simply discarded. The bias introduced by discarding these negative events is easy to evaluate. The absolute value of the discarded negative weight events is superimposed on Figures 6-9 as a shaded histogram. These histograms are 5 orders of magnitude below the 2nd stage integration distributions (too small to resolve on Figures 7 and 9) and are largest in the region where the differential cross section is also largest, and so the effect of these negative events is everywhere negligible.

It is possible to use as few as one sample of each type (n-body and (n+1)-body) for the 2nd stage ‘integration’. The effect of decreasing the samples in this way is to increase the frequency of negative weight events—in this case the bias increases to 1%, which represents a maximum ceiling for the bias.

Finally it is noted that an adaptation of the PSS method may eliminate the negative event weights entirely allowing for the unbiased use of a single sample to evaluate the parton shower region cross section. In Ref.  a hybrid  of the subtraction and the one parameter PSS regularization schemes is investigated. For each n-body configuration the PSS cutoff parameter which gives a zero n-body contribution is calculated—in this manner the virtual and real emission cancellations are enforced and all of the phase space is partitioned to be (n+1)-body. One concern with this method is that the minimum jet energy is coupled to the n-body momentum configuration. However, there can be no bias for any observable which is insensitive to the specific PSS boundary for a given n-body configuration. So long as the hybrid PSS region is everywhere contained within the parton shower region (which is expected to be the case), the average cross section of the parton shower region is one such number—and so this method could be employed to evaluate the 2nd stage integration. In this manner the NLO cross section across the parton shower region (Eq. 5) can be represented with a single sample, which with this hybrid method is always a positive number. The weights of many such events will converge to the correct cross section, and event generation would proceed with a single sampling of the phase space, as it normally does for LO showering event generators. The important difference is that in the parton shower region information about the parton emission would be discarded and the n-body kinematic configuration would be used as input for the parton shower, which eliminates the possibility of the jet energy coupling since, as described previously, the parton shower (guided by the soft matrix element corrections) would determine which exclusive event structure is chosen for the event.

## IV CONCLUSIONS

Matrix element corrections for hadronic diboson production have been implemented following the procedure outlined in Ref.  and are interfaced to the HERWIG showering event generator. The corrections are significant in the high-$`P_{}^T`$ regions where the agreement in shape with NLO calculations is significantly improved.

A new method for incorporating the full NLO matrix element is presented. The inclusive recoil against emissions is dominated by the parton shower for events in the soft and collinear regions where the parton shower is expected to give a good description of the physics, and is dominated by the first order matrix element for events in the hard central region where perturbation theory provides a reliable prediction. The normalization is everywhere NLO. Events generated with this method have positive definite event weight, such that unweighted events can be easily produced by means of Monte Carlo hit-and-miss. The method is generalizable to any color singlet production process.

The numerical implementation of the method is straight-forward providing a NLO matrix element (employing the PSS regularization scheme) is available.

The difference in distribution shapes between the matrix element correction approach and the 2nd stage integration incorporating the full NLO matrix element is small—though noticeable for example in distributions such as the diboson rapidity separation. The 2nd stage integration approach has the further advantage of providing NLO normalization.

Diboson event samples from hadron colliders to date are small and insensitive to QCD corrections, so a meaningful comparison to data will have to await Tevatron Run II and LHC.

## ACKNOWLEDGMENTS

The author would like to thank the ATLAS Collaboration and in particular I. Hinchliffe, M. Lefebvre, and J. B. Hansen. I am grateful to U. Baur for providing the matrix elements, encouragement, and informative correspondence. I benefited from discussions with G. Corcella, M. Mangano, G. Ridolfi, M. Seymour, and B. Webber. I thank L. Bourhis, A. Signer, and J. Stirling for organizing and hosting the 2000 IPPP Workshop on Matrix Elements and Parton Showers in Durham. This work has been supported by the Natural Sciences and Engineering Research Council of Canada.


text:
# Multiply quantized vortices in trapped Bose-Einstein condensates

## I Introduction

Trapped Bose-Einstein condensed gases provide a novel kind of condensed matter system which can sustain quantized vortices, as has been realized in the recent years jila ; ens1 ; ens2 ; mit . Features which differ from those in macroscopic superfluid systems are expected here. The fact that these systems are finite, and that the parameter regime where the vortex core size is comparable to the system size is attainable, in principle opens up for the possibility for vortices in these systems to have a circulation quantum number $`q`$ larger than unity.

It is well known that multiply quantized vortices are not thermodynamically stable in spinless, macroscopic and homogeneous superfluids, because the energy of a vortex depends on the square of the circulation, and therefore two singly quantized vortices with a finite spatial separation have lower energy than one doubly quantized, while giving the system the same angular momentum donnelly . However, multiply quantized vortices may well be energetically favorable in systems which do not fulfill the criteria of being homogeneous, spinless and macroscopic. The assumption of homogeneity is not met in superconductors with pinning forces, where indeed doubly quantized vortices are observed scpinning . In 3He-A, where the order parameter is not a scalar due to spin degrees of freedom, a lattice of doubly quantized vortices with filled cores has recently been observed helsinki . The argument also fails in spatially confined systems where the vortex cores are not much smaller than the system, such as mesoscopic superconducting disks, where vortices with large quantum numbers have been predicted schweigert .

In condensed Bose gases confined in harmonic-oscillator potentials, it has been found both analytically pethick , numerically butts ; castindum  and variationally castindum  that multiply quantized vortices are not energetically favorable, and no quantum numbers larger than unity have been observed in experiment ens2 ; mit . However, none of these studies has considered whether altering the power law of the confining potential can open up for the existence of multiply quantized vortices. In this paper, we shall give rigorous criteria for the thermodynamic stability of multiply quantized vortices in trapped Bose-Einstein condensates, and show the decisive role played by the shape of the potential.

The system under study is a gas of bosons of mass $`m`$ in an external potential $`V(r)`$, dilute enough and at sufficiently low temperature that it is well described by the Gross-Pitaevskii equation gross ; pitaevskii ; thebook

$$i\mathrm{}\frac{\psi (r,t)}{t}=\frac{\mathrm{}^2}{2m}^2\psi (r,t)+V(r)\psi (r,t)+U_0|\psi (r,t)|^2\psi (r,t)\mathrm{\Omega }\widehat{L_z}\psi (r,t).$$
(1)
The so-called condensate wave function $`\psi (r,t)`$ is normalized to the number of particles $`N`$ in three dimensions, and to the number of particles $`\nu `$ per unit length in two dimensions. The coefficient $`U_0`$ in front of the nonlinear term is the interaction strength, defined as $`U_0=4\pi \mathrm{}^2a/m`$, where $`a`$ is the s-wave scattering length. It turns out that the effective measure of interaction strength in two dimensions is the product $`4\pi \nu a`$, which we shall denote by $`g`$. In three dimensions the corresponding quantity depends on the trapping potential.

A centrifugal term is present in the Gross-Pitaevskii equation, corresponding to a rotation of the trap with the frequency $`\mathrm{\Omega }`$ about the $`z`$ axis. At certain critical values $`\mathrm{\Omega }_{\text{c}q}`$ of the rotation frequency, we expect there to be a discrete transition between states of different circulation numbers, so that when $`\mathrm{\Omega }_{\text{c}q+1}>\mathrm{\Omega }>\mathrm{\Omega }_{\text{c}q}`$, a state with total circulation quantum number $`q`$ is the ground state donnelly . The critical frequencies are functions of the coupling strength $`U_0`$ lps .

When the Gross-Pitaevskii equation is valid, the total energy of the system is given by the mean-field energy functional

$$E=𝑑r\psi ^{}(r)\left(\frac{\mathrm{}^2}{2m}^2+V(r)\mathrm{\Omega }\widehat{L}_z\right)\psi (r)+\frac{1}{2}U_0|\psi (r)|^4.$$
(2)
We shall in the following be concerned with the lowest-energy solutions of the Gross-Pitaevskii equation (1) for a range of values of the driving frequency $`\mathrm{\Omega }`$, interaction strength $`U_0`$ and for different external potentials $`V`$. The paper is organized as follows. In Section II, we present an analytic study of systems contained in three- and two-dimensional power-law and anharmonic traps. In Sec. III, we illustrate these findings numerically for the two-dimensional case and study some features of the different states with small and intermediate circulation quantum numbers. Sec. IV provides a conclusion.

## II Vortex configurations in power-law and anharmonic traps

Consider a Bose gas trapped in a cylindrically symmetric power-law potential of power $`n`$, which in cylindrical coordinates is written

$$V(r,\varphi ,z)=\mathrm{}\omega \left(\frac{r}{d_t}\right)^n+V_{}(z).$$
(3)
The trap frequency $`\omega `$, which determines the strength of the potential, defines a trap length $`d_t=\sqrt{\mathrm{}/m\omega }`$. The dependence of the potential on the $`z`$ coordinate is totally arbitrary as long as the problem remains separable in the noninteracting limit, as is the case here. We shall hereafter neglect to mention the irrelevant $`z`$-dependent term and its associated degrees of freedom.

In the analysis of this section, we shall consider the weakly-interacting limit where $`U_0`$ is small, because if a multiply quantized vortex can ever be present, it will be so when the interaction is weak. To see why, note that the vortex cores are large for weak interactions and vice versa thebook ; bp . Therefore, in the limit of large $`U_0`$, the vortices are much smaller than the system size and the radius of curvature of the density profile. The result of a homogeneous bulk system thus applies, namely that a vortex array is the energetically favorable configuration. Therefore, we expect multiply quantized vortices to show up only in the realm of small (and possibly intermediate) $`U_0`$.

Ref. pethick  showed why multiply quantized vortices are never the minimum-energy state in a system which is harmonically confined in the $`x`$-$`y`$ plane. We will here extend the same analysis to the case when the confinement is steeper than harmonic. The argument is based on perturbative treatment of the interaction term in the Gross-Pitaevskii equation (1). We therefore first study the noninteracting gas. When $`U_0=0`$, the Gross-Pitaevskii equation is identical to the one-particle Schrödinger equation, whose solutions are well known when the radial confinement is harmonic. The eigenstates $`\phi _{n_rq}`$ can be labeled by a radial and an axial quantum number $`n_r`$ and $`q`$ (as well as a quantum number associated with the $`z`$ direction, irrelevant for our considerations). The eigenenergies including the centrifugal term are

$$E_{n_rq}=\mathrm{}(\omega \mathrm{\Omega })q+\mathrm{}\omega (1+n_r).$$
(4)
We shall concentrate on the states with no radial nodes, $`\phi _{0q}`$, because these have the lowest energy for a given angular momentum $`Nq\mathrm{}`$. These are in fact vortices with quantum number $`q`$. When the driving frequency $`\mathrm{\Omega }`$ is less than the trap frequency $`\omega `$, the $`q=0`$ state is the ground state; when $`\mathrm{\Omega }=\omega `$, all the states $`\phi _{0q}`$ are degenerate. When $`\mathrm{\Omega }>\omega `$, $`E`$ decreases with increasing $`q`$, and there is no ground state. Hence, in the noninteracting case all the critical frequencies $`\mathrm{\Omega }_{\text{c}q}`$ are equal to the trap frequency $`\omega `$.

Following Ref. pethick , we now consider perturbatively the interaction energy. The perturbation for a given wave function $`\psi `$ is, from Eq. (2),

$$H^{}=\frac{1}{2}U_0𝑑r|\psi (r)|^4.$$
(5)
Calculating this for general superpositions $`\psi (r)=_qC_q\phi _{0q}(r)`$, we find that the degeneracy at $`\mathrm{\Omega }=\omega `$ is lifted, and the critical frequencies $`\mathrm{\Omega }_{\text{c}q}`$ split apart and assume values less than $`\omega `$. Fixing a total angular momentum $`L_z=Nq\mathrm{}`$, it turns out that wave functions which are superpositions of states of different $`q`$ have less interaction energy than the pure multiply-quantized vortex configurations $`\phi _{0q}`$ for $`q2`$. Such superpositions are vortex arrays, and these are the energetically favorable states for any finite interaction strength.

We shall now show how the analysis of Ref. pethick  works when the trapping power $`n`$ in Eq. (3) is larger than two. In this case, the single-particle energy varies faster than linear with the angular momentum. This implies that in the noninteracting case, pure multiply-quantized vortex states $`\phi _{0q}`$ have lower energy than any superposition having the same angular momentum. Hence, it takes a finite (albeit possibly very small) interaction strength $`U_0`$ to overcome this difference. We stress that in order to arrive at this result, we did not need to know any details of the eigenfunctions of a general power-law potential: the fact that the energy depends stronger than linearly on the angular momentum is enough. For a finite but weak enough interaction, the multiply-quantized vortex states are the minimum-energy configurations within their respective ranges of external rotation frequency $`\mathrm{\Omega }`$. When the interaction is strong enough, however, the larger interaction energy of the multiply-quantized vortex states will overcome the difference in kinetic plus trap energy, and the ground state will again be a vortex array. This makes sense, since as we noted above, in the limit of strong interaction the results of the infinite and homogeneous case apply.

For potentials weaker than harmonic, $`n<2`$, no rotation is possible. The potential in a frame rotating with angular frequency $`\mathrm{\Omega }`$ is $`V(r)\frac{1}{2}\mathrm{\Omega }r^2`$, which is not confining if $`\mathrm{\Omega }>0`$ and $`V(r)`$ is weaker than harmonic. Therefore a rotating state in such a trap can be at most metastable. privatePethick

If a pure power-law potential which is steeper than harmonic may seem an artificial construction which is hard to fabricate in practice, we note that the analysis can be repeated for an anharmonic potential,

$$V(r,\varphi ,z)=\frac{1}{2}m\omega ^2r^2(1+\lambda \frac{r^2}{d_t^2})+V_{}(z).$$
(6)
The important point is that for any positive anharmonicity $`\lambda `$, the single-particle energy varies faster than linear with the angular momentum. This again leads to the conclusion that multiply quantized vortices can exist in anharmonic traps as long as $`\lambda >0`$ and the coupling is small.

We have shown that multiply quantized vortices can be stable only in trapping potentials which are steeper than harmonic in the $`x`$-$`y`$ plane, and only for weak enough interactions. A close inspection of the above argument reveals the physical mechanism behind the transition. The kinetic energy gives rise to a repulsion between the vortices because of the overlapping velocity fields, while the trapping potential strives to keep the vortices together in order to reduce the spatial extent of the cloud. The interaction energy favors states with several separated vortices because such states have lower density on the average. The thermodynamic stability of multiply quantized vortices is thus a consequence of the competition between these three energies, and the power 2 is the limiting value.

## III Numerical results for two-dimensional systems

We now perform a numerical study to quantify the predictions of Section II. We shall treat the two-dimensional case only. We insert the power-law potential, Eq. (3), into the Gross-Pitaevskii equation (1) and scale out the dimensions by defining $`r=d_t\stackrel{~}{r}`$, $`t=\omega ^1\stackrel{~}{t}`$, $`\psi =\sqrt{\nu }d_t^1\stackrel{~}{\psi }`$, $`L_z=\mathrm{}\stackrel{~}{L}_z`$, $`\mathrm{\Omega }=\omega \stackrel{~}{\mathrm{\Omega }}`$, $`g=4\pi \nu a`$, and for the total mean-field energy $`E=\mathrm{}\omega \stackrel{~}{E}`$. There results a dimensionless Gross-Pitaevskii equation,

$$i\frac{\stackrel{~}{\psi }}{\stackrel{~}{t}}=\frac{1}{2}_{\stackrel{~}{r}}^2\stackrel{~}{\psi }+\frac{1}{2}\stackrel{~}{r}^n\stackrel{~}{\psi }+g|\stackrel{~}{\psi }|^2\stackrel{~}{\psi }\stackrel{~}{\mathrm{\Omega }}\widehat{\stackrel{~}{L}_z}\stackrel{~}{\psi },$$
(7)
which depends on three parameters, namely the coupling $`g`$, the rotation frequency $`\stackrel{~}{\mathrm{\Omega }}`$, and the trapping power $`n`$. We solve it in two dimensions using the split Fourier method castindum . By propagating Eq. (7) in complex time, the system relaxes to a local energy minimum which depends on the initial state. In order to find the global minimum-energy configuration for each choice of parameters $`g`$, $`\stackrel{~}{\mathrm{\Omega }}`$, $`n`$, we perform the relaxation a few times using different, irregular initial states, compute their energy with the aid of Eq. (2), and out of the different final states we choose that which has the lowest energy. The grid size is $`64\times 64`$ points, and the lattice constant is equal to 0.2 times the trap length $`d_t`$. We have confirmed that doubling the grid size and halving the lattice constant does not alter the results.

Figure 1 is a phase diagram for the quartic potential ($`n=4`$), which shows how the $`\stackrel{~}{\mathrm{\Omega }}`$-$`g`$ plane is divided into regions of different vortex configurations with the total circulation $`q`$ ranging from 1 up to 4. The different vortex configurations are labelled by the total number of quanta followed by the number of singularities present in the system: the configuration 3-3 thus denotes three $`q=1`$ vortices whereas 3-1 is a state with one triply quantized vortex present. The region marked by 0 is the nonrotating state.

Since the numerical computation is done with limited precision, we need an operational criterion for the existence of a multiply quantized vortex. A $`q`$-fold quantized vortex is defined as a configuration in which the calculated circulation around a circle enclosing the origin is equal to $`q`$ and only one density minimum can be seen on optical inspection of the plotted wave function, i. e. the inter-vortex distance is less than the lattice constant of the numerical grid. The transition frequencies found using this criterion coincides with the point where the angular momentum saturates at an integer value.

As expected from the analysis of the previous section, we find multiply quantized vortices for small values of $`g`$, which split into arrays of singly quantized vortices for larger $`g`$. The transition from a multiply quantized vortex to an array is continuous, as illustrated in Figure 2, where two vortices are seen to merge as the critical line is crossed.

It is convenient to define a new quantity to describe the regime of stability of multiply quantized vortices. The $`q`$-th critical coupling $`g_{\text{c}q}`$ is defined so that when $`g>g_{\text{c}q}`$, there can exist no $`q`$-fold quantized vortices for any value of $`\mathrm{\Omega }`$, but when $`g<g_{\text{c}q}`$, a $`q`$-fold quantized vortex is the energy minimum within certain limits of $`\mathrm{\Omega }`$. The critical coupling is a function of the trap power $`n`$. Fig. 1 suggests that the critical coupling increases with $`q`$. The physical reason for this is, that if the system is to accommodate many singly quantized vortices, the cores must be small and $`g`$ large. Our analysis in the previous section showed that $`g_{\text{c}q}=0`$ for the harmonic trap, $`n=2`$, and $`g_{\text{c}q}>0`$ when $`n>2`$. This is confirmed in Figure 3, which shows how $`g_{\text{c}2}`$ depends on $`n`$.

The data points suggest that the curve $`g_{\text{c}2}(n)`$ goes to zero linearly when $`n2`$ from above.

We have pursued the analogous anharmonic case, Eq. (6), using the same numerical procedure. The qualitative features of the phase diagram are not different from the pure power-law case. The critical coupling $`g_{\text{c}q}`$ is now a function of the anharmonicity $`\lambda `$. The result for $`g_{\text{c}2}`$ is shown in Figure 4.

The limit $`\lambda \mathrm{}`$ is the quartic potential, whose critical coupling $`g_{\text{c}2}=30.5`$ can be read off from Fig. 1. When $`\lambda 0`$, on the other hand, the potential approaches a harmonic one, and $`g_{\text{c}2}`$ approaches zero. In the experimentally realized two-dimensional condensates of Ref. mit2d , $`g`$ varies between 100 and $`10^4`$, the former for a number of particles $`N=10^4`$. Diminishing the number of particles by a factor 10 and introducing an anharmonicity larger than 0.1 would take these systems below $`g_{\text{c}2}`$.

The phase diagram in $`g`$-$`\mathrm{\Omega }`$ space becomes richer for higher quantum numbers $`q`$. For $`q7`$, there are only two phases, a vortex-array phase and a multiply-quantized vortex phase. For higher quantum numbers, however, a new kind of pattern appears, where one multiply quantized vortex is surrounded by singly quantized vortices. Figure 5 shows the different vortex phases for $`8q11`$. The power of the potential is chosen to $`n=2.2`$.

We find a number of distinct phases, which we label by their total circulation followed by the number of singularities just as in Fig. 1. There are three phases containing eight quanta of circulation: a central $`q=8`$ vortex for small $`g`$ (denoted 8-1); an array of eight vortices for large $`g`$ (8-8); and for intermediate $`g`$, a central $`q=2`$ vortex surrounded by six $`q=1`$ vortices (8-7). The transition from 8-1 to 8-7 is continuous: the six $`q=1`$ vortices separate from the central one in a smooth way as the critical line is crossed. The two phases 8-7 and 8-8, on the other hand, have different symmetries and the transition between these two is discontinuous. For the case of nine quanta, there exist the corresponding phases 9-1 and 9-9, and one intermediate phase 9-8 with a central $`q=2`$ vortex and seven $`q=1`$ vortices. A 9-7 configuration with a $`q=3`$ vortex surrounded by six $`q=1`$ vortices seems to be stable in a very narrow region of width $`\mathrm{\Delta }g0.2`$, between 9-1 and 9-8, but with the present numerical precision we have not been able to accurately determine the boundaries of this phase and it is left out of Fig. 5. For the case q=10, on the other hand, one can clearly distinguish four different configurations, and likewise for $`q=11`$.

Analogously to the $`q=8`$ case, the transition between 9-1 and 9-7 is continuous, and so is the transition between 10-1 and 10-8 and between 11-1 and 11-8. On the other hand, the configurations 9-7, 10-8 and 11-8 have entirely different symmetries from 9-8, 10-9 and 11-9, respectively, and those transitions are discontinuous. Finally, the three last-mentioned states transform smoothly into arrays of singly quantized vortices (9-9, 10-10 and 11-11) when $`g`$ increases, by splitting of the central vortex into singly quantized ones. Figure 6 illustrates this.

## IV Conclusions

We have determined the criteria for existence of vortices with circulation quantum numbers larger than unity in rotated trapped Bose-Einstein condensates. We have proven that in both two- and three-dimensional systems, multiply quantized vortices are the energy minimum if the trapping potential in the plane perpendicular to the axis of rotation is steeper than harmonic, and the interaction is sufficiently weak so that the vortex cores are not much smaller than the cloud, and the external rotation frequency is within appropriate limits. For stronger interactions, the multiply quantized vortices break up into arrays of several vortices. The physics of the transition is the interplay between the inter-vortex repulsion and the trapping potential. In the case of two-dimensional systems, we have numerically determined the regimes of thermodynamic stability of different vortex configurations. For clouds rotated at large angular velocities, there are regimes where arrays containing both singly and multiply quantized vortices are stable.

Most trapping potentials used in experiments are harmonic, but using Laguerre-Gaussian beams one can construct optical traps of arbitrary even power laws, where our predictions can be tested kuga . A two-dimensional Sodium condensate can host doubly quantized vortices if the number of particles is of the order 1000 or smaller, and a trap anharmonicity of about 10% is introduced. Alternatively, the appropriate weak-coupling regime can be attained by using Feshbach resonances to diminish the scattering length feshbach . In a three-dimensional system, a quantitative determinition of the conditions has been prohibited by numerical limitations. The technique of vortex imaging by interference interference , which directly images the phase of the condensate, promises to be a convenient way to detect multiply quantized vortices. Alternatively, one can look for multiply quantized vortices by measuring the angular momentum of the system angmom , or simply by looking for anomalously large vortex cores.

###### Acknowledgements.

I am grateful to my supervisor Jørgen Rammer for supervising this work, and Ping Ao for supervision at the early stages. I would also like to thank Georgios Kavoulakis, C. J. Pethick, Thomas Busch and Lars Melwyn Jensen for valuable discussions, NORDITA for travel support and lodging, and Jani Martikainen for help with the numerics.