text:
# Mesoscopic superconductors in the London limit: equilibrium properties and metastability

## 1 Introduction

A significant amount of work has been recently devoted to the study of aluminum superconducting disks in a mesoscopic regime, in which the size of the sample is comparable with both the coherence length $`\xi `$ and the London penetration depth $`\lambda `$. In a first set of experiments  the magnetization of such systems was measured at thermodynamic equilibrium. Beyond the Meissner state, a series of discontinuous jumps appear when the applied magnetic field is increased. This behaviour corresponds to the entrance of individual quantized vortices in the sample. A quantitative understanding of the phenomena involved can be achieved through a numerical study of the Ginzburg-Landau equations . Analytical progress can be made if one considers special limits of these equations (such as the linear limit , the London regime  or the Bogomol’nyi dual point ). In later experiments , it was noticed that when sweeping down the applied field the sample exhibits a paramagnetic Meissner effect, i.e. it has a paramagnetic magnetization. Such a behaviour previously found in high $`T_c`$ superconductors  had been interpreted as a special features of these materials. However, its occurrence in aluminum disks calls for a less specific explanation which results from the role played by metastable states at temperatures well below the critical temperature.

This work is devoted to the study of metastable states in a mesoscopic superconductor. We shall investigate vortex patterns in the London regime, i.e. in the limit $`\kappa 1`$, where $`\kappa `$ is the Ginzburg-Landau parameter. Using known results  for the magnetization of an infinite cylinder with a circular cross-section of radius $`R`$, we derive, in the mesoscopic regime, closed expressions for the magnetic field and the free energy that are suitable for analytical calculations. They allow to obtain the Bean-Livingston barrier and to calculate the series of matching fields (critical fields corresponding to the entrance of the $`N`$-th vortex). The paramagnetic Meissner state can then be explained using a simple scenario for the hysteretic behaviour of the magnetization, involving the appearance and the disappearance of the surface barrier.

Once the existence of metastable states established, a criterion to classify them is required. The number of vortices is a simple topological number but it is not precise enough to distinguish between all different type of states. Projection methods into eigenstates of the linearized theory produce a set of integers that were proposed as classifying numbers to characterize a given metastable vortex state  and that would generalize quantum numbers to a non-linear theory. However, these numbers are non-generic and it is not obvious to relate them to the geometric features of a given pattern. A more geometrical path will be followed here to present a topological description of vortex patterns. Our study is based on an analogy between a configuration of vortices and a dynamical system, obtained by interpreting the superconducting current as the local (phase-space) velocity of a particle. The Hamiltonian of the associated dynamical system turns out to be the magnetic field. We shall show that a configuration of vortices can be characterized by the critical points of the magnetic field (or equivalently the points where the current vanishes). The number of critical points of a given type (maxima, minima, saddle-points) will provide a natural set of topological invariant integers associated with a configuration. This geometrical analysis reveals the existence of structural transitions between states with the same number of vortices: indeed, the number and the nature of the critical points can vary abruptly when the applied field exceeds some specific values, provoking topological phase transitions. For a simple polygonal configuration of vortices, we calculate the corresponding transition fields. These topological changes are in fact best visualized using a ‘dual’ description of a system of vortices, closely related to the paramagnetic Meissner effect: there exists a special closed curve $`\mathrm{\Gamma }`$, first introduced in , which is everywhere orthogonal to the current lines and which passes through the critical points of the system. Therefore any structural change affecting the critical points, will manifest itself in a topological change of the shape of $`\mathrm{\Gamma }`$. This curve thus provides an efficient and elegant tool to characterize vortices patterns and to understand qualitatively the topological transitions between different configurations.

The plan of this article goes as follows. In section 2, we recall how the London free energy can be obtained as a limit of the Ginzburg-Landau theory. In section 3, we present the solution of London equation in the case of an infinite circular cylinder (details of the calculations are given in Appendix A and B). Section 4 is devoted to the mesoscopic limit: we calculate the magnetic field and the free energy, study the Bean-Livingston barrier and obtain the matching fields. The paramagnetic Meissner state is obtained in section 5. Section 6 is devoted to the topological investigation of vortex patterns: we show that the critical points characterize a configuration and we demonstrate the existence of topological phase transitions. In the second part of section 6, the dual description of a configuration of vortices, via the curve $`\mathrm{\Gamma },`$ is introduced to provide a better understanding of the transitions between patterns. Section 7 contains some concluding remarks and discussions.

## 2 The London limit of the Ginzburg-Landau equations

Consider an infinite superconducting cylinder with a cross-section of area $`\mathrm{\Omega }`$ in the constant magnetic field $`h_e\widehat{z}`$ along the axis of the cylinder. In the presence of a magnetic field the Ginzburg-Landau free energy density per unit length of the cylinder, defined as difference of free energies with and without magnetic field $`=F_S(h)F_S(0)`$, reads

$$=\left(\frac{\stackrel{}{h}^2}{2}+\kappa ^2\left(1|\psi |^2\right)^2+\left|\left(\stackrel{}{}i\stackrel{}{A}\right)\psi \right|^2\right)d^2r$$
(1)
where $`\psi =|\psi |e^{i\chi }`$ is the order parameter and $`\stackrel{}{h}=\stackrel{}{}\times \stackrel{}{A}`$ is the local magnetic induction. The integration is performed over the cross-section of the cylinder. A superconductor is characterized by two length scales: the London penetration depth $`\lambda `$ and the coherence length $`\xi `$. The Ginzburg-Landau parameter is defined as their ratio $`\kappa =\lambda /\xi `$. In this work lengths are measured in units of $`\lambda \sqrt{2}`$ and the magnetic induction in units of the thermodynamic critical field $`H_c=\varphi _0/2\sqrt{2}\pi \lambda \xi `$, where $`\varphi _0=hc/2e`$ is the magnetic flux quantum. The free energy (1) is given in units of $`\xi ^2H_c^2/4\pi `$.

The Ginzburg-Landau equations for $`\psi `$ and $`\stackrel{}{h}`$ are obtained from a variation of $``$. They are given by

$`\left(\stackrel{}{}i\stackrel{}{A}\right)^2\psi `$ $`=`$ $`2\kappa ^2\psi \left(1|\psi |^2\right)`$ (2)
$`\stackrel{}{}\times \stackrel{}{h}`$ $`=`$ $`2|\psi |^2\left(\stackrel{}{}\chi \stackrel{}{A}\right)`$ (3)
Equation (3) is the Maxwell-Ampère equation which could be written as well using the current density $`\stackrel{}{ȷ}=\text{Im}\left(\psi ^{}\stackrel{}{}\psi \right)|\psi |^2\stackrel{}{A}`$. Outside the superconducting sample, $`\psi =0.`$ The boundary condition on the surface of the superconductor is obtained by requiring that the normal component of the current density vanishes (superconductor/insulator boundary condition ):

$$\left(\stackrel{}{}i\stackrel{}{A}\right)\psi |_{\widehat{n}}=0$$
(4)
here $`\widehat{n}`$ is the unit vector normal at each point to the surface of the superconductor.

We study a superconductor in an applied external field $`h_e`$ the relevant thermodynamic potential is the Gibbs free energy $`G`$ obtained from $`F`$ via a Legendre transformation

$$G=F\mathrm{\Omega }h_eB$$
(5)
where the total (dimensionless) magnetic induction (in the $`\widehat{z}`$ direction ) $`B`$ is given by the averging the magnetic induction over the cross-section

$$B=\frac{1}{\mathrm{\Omega }}h(\stackrel{}{r})d^2r$$
(6)
where the magnetic induction is always in $`\widehat{z}`$ direction: $`\stackrel{}{h}=h\widehat{z}`$.

In the normal phase, we have $`\psi =0`$, $`B=h_e`$ and $`F_S(h)=F_N(h)=\mathrm{\Omega }h^2/2`$. Therefore, the corresponding Gibbs free energy $`G_N`$ is given by:

$$G_N=F_N(h)\mathrm{\Omega }h_eB=F_N(0)\frac{\mathrm{\Omega }h_e^2}{2}$$
(7)
At thermodynamic equilibrium, the superconductor selects the state of minimal Gibbs free energy $`𝒢`$. The quantity which is measured in experiments is the magnetization $`M`$ of the superconductor due to the applied field given by $`4\pi M=Bh_e`$. It is obtained, at thermodynamic equilibrium using the thermodynamic relation :

$$M=\frac{1}{2\pi }\frac{𝒢}{h_e}$$
(8)
where the difference of the (dimensionless) Gibbs energies $`𝒢=G_SG_N`$ up to a constant equal to the superconducting condensation energy, is given by

$$𝒢=\mathrm{\Omega }h_eB+\mathrm{\Omega }h_e^2/2$$
(9)
The dimensionless ratio $`\kappa =\lambda /\xi `$ is the only free parameter to describe the superconducting state. It determines, in the limit of an infinite system, whether the sample is a type-I or type-II superconductor . For $`\kappa 1/\sqrt{2}`$, i.e. for type-II superconductors, it is energetically favorable for the system in a large enough magnetic field to sustain normal regions in the bulk which appear as vortices. Their magnetic flux quantized in units of $`\varphi _0`$. We shall now study the system of vortices in an extreme Type II superconductor, i.e. in the London limit $`\kappa 1`$.

In this limit, the amplitude of the order parameter $`|\psi |`$ does not vary and can be taken to one except at the position of the vortices where it vanishes. The Ginzburg-Landau free energy reduces to

$$=\left(\frac{h^2}{2}+\left|\left(\stackrel{}{}i\stackrel{}{A}\right)\psi \right|^2\right)d^2r$$
(10)
Using this equation and $`|\psi |^2=1`$, except at the position of the vortices, we can write

$$\left|\left(\stackrel{}{}i\stackrel{}{A}\right)\psi \right|^2=\left(\stackrel{}{}\times \stackrel{}{h}\right)^2$$
(11)
so that the free energy is given by the expression

$$=\frac{1}{2}\left(h^2+\frac{1}{2}\left(\stackrel{}{}\times \stackrel{}{h}\right)^2\right)d^2r$$
(12)
which in our units is independent of $`\kappa `$. To obtain the equation for $`h`$, we notice that the phase $`\chi `$ is non singular except near a vortex (say at $`r=0`$) where it has the property that

$$\stackrel{}{}\chi \stackrel{}{d}l=2\pi $$
(13)
on any closed contour encircling the singularity at point $`r=0`$. Using Stokes theorem,

$$\stackrel{}{}\chi \stackrel{}{d}l=𝑑\stackrel{}{S}\stackrel{}{}\times \stackrel{}{}\chi =2\pi ,$$
(14)
we deduce that the $`\widehat{z}`$-component of $`\stackrel{}{}\times \stackrel{}{}\chi `$ equals $`2\pi \delta (\stackrel{}{r})`$. The equation for $`h`$ is then given by

$$^2h2h=4\pi \delta (\stackrel{}{r})$$
(15)
This equation can be generalized to the case of $`N`$ vortices placed at points $`\stackrel{}{r}_k`$

$$^2h2h=4\pi \underset{k=1}{\overset{N}{}}\delta (\stackrel{}{r}\stackrel{}{r}_k)$$
(16)
In the next section we shall present the solution of this equation in the particular case of a cylinder with circular cross section.

## 3  Solution for $`N`$ vortices in a circular cylinder

¿From now on we shall study a circular cylinder of radius $`R`$. The magnetic induction and Gibbs energy of a system of $`N`$ vortices in this geometry has been studied first in . The detailed derivation of the results is presented in the Appendix A. The solution $`h(r,\theta )`$ of the equation (16) satisfying the boundary condition $`h(R,\theta )=h_e`$ can be written as a sum of three terms:

$$h(r,\theta )=h_M(r)+h_V(r,\theta )+h_{\overline{V}}(r,\theta )$$
(17)
The first term

$$h_M(r)=h_eI_0(\sqrt{2}r)/I_0(\sqrt{2}R),$$
(18)
where $`I_0(x)`$ is the modified Bessel function of first kind, describes the Meissner effect in the absence of vortices. The second term is the magnetic induction of $`N`$ vortices placed at points $`\stackrel{}{r}_k=(r_k,\theta _k)`$ for $`k=1,\mathrm{},N`$ and it is given by

$$h_V(r,\theta )=2\underset{k=1}{\overset{N}{}}K_0(\sqrt{2}|\stackrel{}{r}\stackrel{}{r}_k|)$$
(19)
The third term, written with help of the modified Bessel function of first and second kind $`I_n(x)`$, $`K_n(x)`$ as

$$h_{\overline{V}}(r,\theta )=2\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}K_n(\sqrt{2}R)\frac{I_n(\sqrt{2}r)}{I_n(\sqrt{2}R)}\underset{k=1}{\overset{N}{}}I_n(\sqrt{2}r_k)\mathrm{cos}n(\theta \theta _k)$$
(20)
ensures the boundary condition $`h(R,\theta )=h_e`$. We are interested in polygonal rings of vortices. For this particular configuration $`r_k=d`$, $`\theta _k=2\pi k/N`$ and the expression for $`h_{\overline{V}}(r,\theta )`$ can be simplified

$$h_{\overline{V}}(r,\theta )=2N\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}K_n(\sqrt{2}R)\frac{I_n(\sqrt{2}r)I_n(\sqrt{2}d)}{I_n(\sqrt{2}R)}\mathrm{cos}Nn\theta $$
(21)
The corresponding Gibbs energy has also been obtained in . The expression of the dimensionless Gibbs energy is given by

$$𝒢=h_eh(\stackrel{}{r})d^2r+\pi R^2h_e^2/2$$
(22)
The field $`h_V`$ diverges at the positions of the vortices and the corresponding part of the free energy must be regularized. The singular contribution coincides with one calculated by Abrikosov  for the case of infinite system and is given by

$$_{\mathrm{}}=N+2\pi \underset{jk}{}K_0(\sqrt{2}|\stackrel{}{r}_j\stackrel{}{r}_k|)$$
(23)
where $`=2\pi (\mathrm{ln}\kappa +0.081)`$ is the one-vortex energy.

The regular part of the Gibbs energy is calculated in the Appendix B. Using the expressions (18),(19) and (20) for the magnetic field, one obtains

$`𝒢_N𝒢_0`$ $`=`$ $`N+2\pi {\displaystyle \underset{jk}{}}K_0(\sqrt{2}|\stackrel{}{r}_j\stackrel{}{r}_k|)2\pi h_e{\displaystyle \underset{k=1}{\overset{N}{}}}\left(1{\displaystyle \frac{I_0(\sqrt{2}r_k)}{I_0(\sqrt{2}R)}}\right)`$ (24)
$``$ $`2\pi {\displaystyle \underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}{\displaystyle \frac{K_n(\sqrt{2}R)I_n(\sqrt{2}r_i)I_n((\sqrt{2}r_j)}{I_n(\sqrt{2}R)}}\mathrm{cos}n(\theta _i\theta _j)`$
where $`𝒢_0=\pi Rh_e^2\left(\frac{R}{2}\frac{\sqrt{2}}{2}\frac{I_1(\sqrt{2}R)}{I_0(\sqrt{2}R)}\right)`$ is the Gibbs energy in the absence of vortices i.e. in the Meissner state. For the particular case of vortices distributed on a regular polygon of radius $`d`$, the previous expression (24) simplifies and becomes

$`{\displaystyle \frac{𝒢_N𝒢_0}{N}}`$ $`=`$ $`+2\pi {\displaystyle \underset{l=1}{\overset{N1}{}}}K_0(|2\sqrt{2}d\mathrm{sin}{\displaystyle \frac{l\pi }{N}}|)2\pi h_e\left(1{\displaystyle \frac{I_0(\sqrt{2}d)}{I_0(\sqrt{2}R)}}\right)`$ (25)
$``$ $`2\pi N{\displaystyle \underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{K_{Nn}(\sqrt{2}R)I_{Nn}(\sqrt{2}d)^2}{I_{Nn}(\sqrt{2}R)}}`$
Similar expressions have been already obtained previously in  and we shall now discuss them in the limit of a small radius $`R`$ in order to use the particularly elegant structure of the corresponding solutions to derive analytical expressions for the field and the energy.

## 4 The limit of a small disk

We shall now consider the limiting case where $`R\lambda `$ for which the corresponding Gibbs free energy was obtained in . In this limit the previous expressions greatly simplify. Since all the distances are given in units of $`\lambda \sqrt{2}`$, we can expand the expression for the magnetic induction using the expansion of the Bessel functions for small argument, so that

$$h(z,\overline{z})h_e=\varphi _e(|z|^21)+2\underset{i=1}{\overset{N}{}}\mathrm{ln}\left|\frac{1\overline{z_i}z}{zz_i}\right|$$
(26)
where $`\varphi _e=h_eR^2/2`$ is the flux of the external field through the system and we use the complex notations,

$$z=\left(\frac{r}{R}\right)e^{i\theta }\text{ and }z_k=\left(\frac{r_k}{R}\right)e^{i\theta _k}$$
(27)
together with complex conjugate $`\overline{z}`$ and $`\overline{z}_k`$. This expression for the magnetic induction is solution of the modified London equation

$$\frac{^2}{z\overline{z}}h(z,\overline{z})h_e=\pi \underset{i=1}{\overset{N}{}}\delta (zz_i)\delta (\overline{z}\overline{z}_i)$$
(28)
which can be obtained from (16) for small variation of the field on the scale of $`R`$.

In the particular case, we shall consider from now on, of vortices distributed along a regular polygon of radius $`d`$, we have $`r_k=d`$ and for the $`N`$ vortices sitting on a shell, the angle $`\theta _k`$ of the $`k`$-th vortex is $`\theta _k=2\pi k/N`$. Then, defining the real quantity $`x=d/R1`$, the relation (26) becomes

$$h_N(z,\overline{z})h_e=\varphi _e(|z|^21)+2\mathrm{ln}\left|\frac{1x^Nz^N}{x^Nz^N}\right|$$
(29)
whereas the Gibbs energy (25) becomes

$$\frac{𝒢_N𝒢_0}{N}=^{}2\pi \varphi _e(1x^2)+2\pi \mathrm{ln}(1x^{2N})2\pi (N1)\mathrm{ln}x2\pi \mathrm{ln}N$$
(30)
where $`𝒢_0=\pi \varphi _e^2/2`$ and $`^{}2\pi \mathrm{ln}\left(R/\xi \right)`$.

The relation between the radius $`x`$ of the polygonal configuration and the external flux $`\varphi _e`$ is obtained by minimizing $`𝒢_N(x,h_e)`$ with respect to $`x`$ at fixed $`\varphi _e`$, namely

$$\varphi _e=\frac{N1}{2x^2}+\frac{Nx^{2N2}}{1x^{2N}}$$
(31)
### 4.1 The Bean-Livingston barrier

The relation (31) has, as a function of the applied flux $`\varphi _e`$, either zero, one or two solutions. The latter case corresponds for the Gibbs energy (30) to the existence of a Bean-Livingston potential barrier  for which one solution is stable and gives the equilibrium position of the vortices while the second solution is unstable and gives the height of the potential barrier as shown in the Fig. 1. It should be noticed that for more than one vortex $`(N>1)`$, the energy $`𝒢_N`$ diverges logarithmically at the origin $`x=0`$ due to the repulsion between the vortices.

The potential barrier disappears for low enough applied flux below a characteristic value $`\varphi _{min}`$ given by the minimum of the function $`f_N(x)`$ in the rhs of (31) as represented in Fig. 2. Defining $`y=x^2`$, we have

$`f_N(y)`$ $`=`$ $`{\displaystyle \frac{N1}{2y}}+{\displaystyle \frac{Ny^{N1}}{1y^N}}`$ (32)
$`f_N^{}(y)`$ $`=`$ $`{\displaystyle \frac{N1}{2y^2}}+{\displaystyle \frac{N(N1)y^{N2}}{1y^N}}+{\displaystyle \frac{N^2y^{2N2}}{(1y^N)^2}}`$ (33)
For large values of $`N`$, we have $`y_{min}^N1/2N`$ so that

$$\varphi _{min}\frac{N1}{2\sqrt[N]{2N}}\frac{N1}{2}$$
(34)
For $`\varphi _e<\varphi _{min}`$, there is no stable solution and then no equilibrium position for $`N`$ vortices exists.

### 4.2 Matching fields

By plugging the expression (31) into the expression (30), we obtain the behavior of the Gibbs energy for $`N`$ vortices as a function of the applied flux $`\varphi _e`$. This expression is the counterpart, in the London limit, of the parabolae obtained at the Bogomoln’yi point  for $`\kappa =1/\sqrt{2}`$. For a given value $`\varphi _e`$ of the applied field, the optimal number of vortices corresponds to the envelop of the set of curves $`𝒢_N(\varphi _e)`$ as represented in the Fig. 3. The matching values $`\varphi _N`$ are the values taken by the flux of the external field at which $`𝒢_{N1}(\varphi _N)=𝒢_N(\varphi _N)`$. They correspond to the transition between configurations with $`N1`$ and $`N`$ vortices respectively.

Using the relation (30), it is straightforward to obtain the field at which the first vortex penetrates the system; it is given by $`\varphi _1=^{}/2\pi =\mathrm{ln}\left(R/\xi \right)`$. This field is larger than the field at which the Bean-Livingston potential barrier disappears. For larger values of $`N`$, in the limit $`\varphi _e\varphi _{min}`$, we approximate the relation (31) by $`\varphi _e=(N1)/2x^2`$ so that

$$\frac{𝒢_N𝒢_0}{2\pi }N\left(\varphi _1\varphi _e\right)+\frac{N(N1)}{2}\left(1+\mathrm{ln}\frac{2\varphi _e}{N1}\right)N\mathrm{ln}N$$
(35)
Equating this to the similar expression for $`𝒢_{N1}(\varphi _e)`$ we obtain that the matching fields for $`N\mathrm{}`$ are given by solution of the following equation

$$\varphi _N=\varphi _1+\frac{N1}{2}\left(1+2\mathrm{ln}\frac{2\varphi _N}{N}\right)\mathrm{ln}N+1+𝒪(1/N)$$
(36)
We use the following ansatz $`\varphi _N=\frac{1}{2}(aN+b\mathrm{ln}N+c)`$. Substituting it to (36) we obtain the coefficients

$$b=\frac{2a}{a2},c=\frac{2a}{a2}\left(1+\varphi _1\frac{a}{2}\right),$$
(37)
where the coefficient $`a`$ satisfies the transcendental equation:

$$a=1+2\mathrm{ln}a$$
(38)
This equation has one obvious solution $`a=1`$ which is compatible with (36). However from the exact numerical calculations of the matching fields using the relation (31) we infer that the second root $`a3.5129`$ is realized in the solution for the matching fields. It leads to the formula

$$\varphi _N=1.76N2.32\mathrm{ln}N+2.32\varphi _11.76+𝒪((\mathrm{ln}N)^2/N)$$
(39)
In the Fig. 4 the values of the matching fields obtained from the computed free energy are compared with our prediction. The discrepancy is attributed to the slow convergence due to the term $`𝒪((\mathrm{ln}N)^2/N)`$.

### 4.3 Magnetization and paramagnetic Meissner effect

A paramagnetic behavior for the total magnetization results from the existence of a high enough Bean-Livingston barrier and the long life-time of the corresponding metastable states. We propose the following scenario for the hysteretic behavior of the magnetization when the magnetic field is increased and then swept down.

If thermal equilibrium is maintained while increasing slowly the external flux $`\varphi _e`$, the results of the preceding section hold and the magnetization is given by the derivative of the envelop of the curves $`𝒢_N(\varphi _e)`$ according to the thermodynamic relation (8). This behavior corresponds to the upper curve in the Fig. 5.

When the applied field is decreased, the vortices that are already in the system are confined by the Bean-Livingston barrier and cannot escape the sample until the flux $`\varphi _e`$ reaches a certain minimum value at which the barrier can be overcome, for instance, by the effect of thermal fluctuations. For the present discussion we use the ultimate criterion of zero barrier as a definition of the transition field.

¿From the discussion of section 4.1, we deduced that the potential barrier for $`N`$ vortices disappears for a characteristic field given for large $`N`$ by the relation (34), $`\varphi _{min}(N1)/2`$. The number of vortices is then decreased by one at this field and the energy drops as shown by the arrows in Fig. 3. The magnetization calculated from (8) is shown by the lower curve in Fig. 5. It shows clearly a paramagnetic behavior, although it disagrees quantitatively with experiments and the numerical simulations based on the exact solution of the Ginzburg-Landau equations. This discrepancy results from the rather unrealistic London limit we used for systems for which $`\kappa 0.3`$ in order to obtain an analytic expression for the Gibbs energy. We emphasize that the system is not at thermal equilibrium but in a metastable state. It is worth noticing that for values of the parameter $`\varphi _1=^{}/2\pi =\mathrm{ln}R/\xi 2`$, the paramagnetic Meissner effect can be observed as well for an equilibrium configuration, i.e. for a minimum of the Gibbs free energy. But this corresponds to $`R\xi `$ so that the delta-like vortex approximation, crucial for the present discussion, is not justified.

## 5 Topological phase transitions between metastable vortex patterns

Up to now, we have considered vortex configurations for which the position $`x`$ of the vortex ring was given by (31) resulting from the minimization of the energy. We shall now relax this condition and study the behavior of the metastable vortex patterns as a function of the applied magnetic flux $`\varphi _e`$. Physically, this corresponds to a situation where $`N`$ vortices are pinned by a potential and remain on fixed positions.

### 5.1 An effective Hamiltonian system

When varying the field $`\varphi _e`$, transitions in a vortex pattern with fixed number $`N`$ of vortices will appear. To characterize them, we shall study the behavior of the magnetic field $`h`$ given by (29). To that purpose, it is interesting first to notice that, in two dimensions, the Maxwell-Ampère equation $`\stackrel{}{ȷ}=\stackrel{}{}\times h\widehat{z}`$ is the Hamilton equation of a system whose generalized coordinates $`(p,q)`$ are the coordinates $`(x,y)`$ such that $`j_x=\dot{x}`$, $`j_y=\dot{y}`$ and the canonically conjugated momentum:

$$\{\begin{array}{ccc}\dot{p}& =& \hfill _qh\\ \dot{q}& =& \hfill _ph\end{array}.$$
(40)
Then, the magnetic field corresponds to the Hamiltonian of the system and the current lines are the phase space trajectories. For a Hamiltonian system, the flow can be characterized by its critical (fixed) points at which

$$\stackrel{}{}h=0$$
(41)
Using (40), we see that critical point correspond to zero velocity points such that $`\stackrel{}{ȷ}=0`$.

In such a description, a vortex corresponds to a maximum in phase-space. However there also exist in the system other critical points, such as minima and saddle-points. The number $`N_M=N`$ of maxima (i.e. vortices), $`N_m`$ of minima, and $`N_s`$ of saddles are not independent, because of the Euler-Poincare relation , namely the Euler-Poincaré characteristic $`\chi `$ is given by

$$\chi =N_MN_s+N_m$$
(42)
The integer $`\chi `$ is a topological invariant and is equal to 1 for a disk thus providing a constraint between the numbers of the different types of critical points. Hence, for a system with a given number of vortices, the difference $`N_sN_m`$ between saddles and minima is fixed. However, each of them can vary and the set of these numbers provides a complete description of the topology of the vortex configurations, namely each topological phase can be described using this set of integers.

A different way to encode the position and the distribution of the critical points of the magnetic field in the disk is to use a special contour $`\mathrm{\Gamma }`$, introduced in . The curve $`\mathrm{\Gamma }`$ is defined by the condition that at each point $`\stackrel{}{r}`$ it is perpendicular to the current density $`\stackrel{}{ȷ}(\stackrel{}{r})`$. The equation of $`\mathrm{\Gamma }`$ is then

$$\frac{dr}{d\theta }=r^2\frac{_rh}{_\theta h}$$
(43)
To this definition we must add the requirement that $`\mathrm{\Gamma }`$ must pass through the critical points at which $`\stackrel{}{ȷ}(\stackrel{}{r})=0`$. The curve $`\mathrm{\Gamma }`$ has several branches; one of them encircles all the vortices and defines a natural boundary between the diamagnetic (Meissner) domain near the boundary and the inner paramagnetic domain which includes the vortices (Fig. 6). An example of this curve is shown in Fig. 7 in the case of a disk with one vortex.

The study of the critical points and the resulting behavior of $`\mathrm{\Gamma }`$ provides a complete description of the vortex states either for thermodynamic or metastable states and of the transitions between different sets of topological numbers when varying the applied field. ¿From the behavior of the curve $`\mathrm{\Gamma }`$ as a function of the applied field $`\varphi _e`$, it might be interesting to draw the analogy with the pressure exerted on a closed membrane separating two systems.

### 5.2 Topological study of a vortex ring

We now consider the simple case where all the vortices are placed on regular polygon at a distance $`x`$ from the center of the disk. The critical points can be classified (see Appendix C) according to their stability properties, obtained from the hessian matrix $`H_{ij}=_i_jh`$. We review the various topological structures obtained and the transitions between them.

For a ring configuration of $`N`$ vortices located at the distance $`x`$ from the center the general expression (26) for the field $`h`$ in circular coordinates $`r`$, $`\theta `$ reads

$$h=\varphi _e(r^21)+\mathrm{ln}\left(\frac{12x^Nr^N\mathrm{cos}N\theta +x^{2N}r^{2N}}{x^{2N}2x^Nr^N\mathrm{cos}N\theta +r^{2N}}\right)$$
(44)
At the critical point the partial derivatives of $`h`$ vanish. The condition $`_rh=0`$ leads to

$$\varphi _erN(x^{2N}1)r^{N1}\frac{x^N\mathrm{cos}N\theta (r^{2N}+1)(x^{2N}+1)r^N}{\left(12x^Nr^N\mathrm{cos}N\theta +x^{2N}r^{2N}\right)\left(x^{2N}2x^Nr^N\mathrm{cos}N\theta +r^{2N}\right)}=0$$
(45)
while the condition $`_\theta h=0`$ rewrites as

$$Nr^N\mathrm{sin}N\theta \frac{(1x^{2N})(1r^{2N})}{\left(12x^Nr^N\mathrm{cos}N\theta +x^{2N}r^{2N}\right)\left(x^{2N}2x^Nr^N\mathrm{cos}N\theta +r^{2N}\right)}=0$$
(46)
The solutions of equations (45) and (46) are studied in detail in Appendix C. One main feature is that one has to distinguish between the cases $`N=1`$, $`N=2`$ and $`N>2`$ where $`N`$ is the number of vortices. We now describe these three cases separately.

#### 5.2.1 Metastable configurations for $`N=1`$

For one single metastable vortex located at distance $`x`$ from the center there is no other critical point for weak flux $`\varphi _e<\varphi _m(1)`$, with

$$\varphi _m(1)=\frac{1x}{1+x}.$$
(47)
At $`\varphi _e=\varphi _m(1)`$ a critical point appears at the point $`r=1,\theta =\pi `$. For a flux slightly above $`\varphi _m(1)`$, this point splits into a minimum at $`r<1,\theta =\pi `$ and two saddle points on the boundary $`r=1`$ situated symmetrically about the diameter passing through the vortex. As the flux increases the minimum moves towards the center and the saddle points approach the point $`r=1,\theta =0`$ each from its side. They coalesce at this point for $`\varphi _e=\varphi _M(1)`$, where

$$\varphi _M(1)=\frac{1+x}{1x}$$
(48)
and for larger values of the flux one saddle point appears inside the disk between the vortex and the boundary in the interval $`x<r<1,\theta =0`$. When $`\varphi _m(1)<\varphi _e<\varphi _M(1)`$, the definition (42) must be generalized to the case of critical point on the boundary . It becomes

$$\chi =N_MN_s+N_m\frac{1}{2}N_{sb}$$
(49)
where $`N_{sb}`$ is the number of saddle points on the boundary. Thus the index is again equal one.

An efficient way to visualize these transitions is examine the changes in the curve $`\mathrm{\Gamma }`$ as the magnetic flux varies. In the case $`\varphi _e>\varphi _M(1),`$ shown in Fig. 7a, the curve $`\mathrm{\Gamma }`$ lies inside the disk and separates the paramagnetic (internal) and diamagnetic (external) regions. For $`\varphi _e=\varphi _M(1)`$ the saddle point reaches the boundary of the disk and for $`\varphi _m(1)<\varphi _e<\varphi _M(1)`$ the curve $`\mathrm{\Gamma }`$ starts and ends on the saddle points on the boundary (see Fig. 7b). Therefore, in this regime, one cannot define a paramagnetic region inside the disk and there is no closed curve that encloses a unit of flux. There appears an arc of circle on the boundary where the direction of the current flow is opposite to the screening current responsible for the Meissner effect. At $`\varphi _e=\varphi _m(1)`$ all the disk becomes paramagnetic and the curve $`\mathrm{\Gamma }`$ disappears.

#### 5.2.2 Metastable configurations for $`N=2`$

For low enough field $`\varphi _e<\varphi _m(2)`$, with

$$\varphi _m(2)=2\frac{1x^2}{1+x^2}.$$
(50)
the only critical point is a saddle point at the center. For $`\varphi _m(2)<\varphi _e<\varphi _M(2)`$, where

$$\varphi _M(2)=2\frac{1+x^2}{1x^2}.$$
(51)
there are two minima situated symmetrically on the diameter perpendicular to the line passing through the vortices and two pairs of saddle points on the boundary. This configuration is shown in Fig. 8a. Each pair is symmetric with respect to the line passing through the vortices and moves towards the points $`r=1,\theta =0`$ or $`\theta =\pi `$. For $`\varphi _e>\varphi _M(2)`$ each pair disappears and gives birth to one single saddle point located between a vortex and the boundary at $`\theta =0`$ or $`\pi `$ as shown in Fig. 8c.

Two minima inside the disk move towards the center and at the critical value of the flux $`\varphi _e=\varphi _c(2)`$, where $`\varphi _c(2)=(1x^4)/x^2`$ they merge with the central point, which becomes a minimum for $`\varphi _e>\varphi _c(2)`$. This critical value can be shown to be always greater than $`\varphi _m(2)`$, but it can be greater or less than $`\varphi _M(2)`$ depending on the value of $`x`$. For $`x<\sqrt{2\sqrt{3}}0.52`$, one has $`\varphi _c(2)>\varphi _M(2)`$, while for $`x>\sqrt{2\sqrt{3}},`$ we have $`\varphi _c(2)<\varphi _M(2)`$. The latter case is shown in Fig. 8b and 8c.

#### 5.2.3 The case $`N>2`$: topological phase transitions

In that case, the central point is a minimum for all values of $`\varphi _e`$. In addition for low enough external flux there are $`N`$ saddle points inside the disk situated symmetrically at $`\theta =2\pi (k+1/2)/N`$. As the external flux is increased through $`\varphi _m(N)`$, where

$$\varphi _m(N)=N\frac{1x^N}{1+x^N}.$$
(52)
$`N`$ pairs of saddle points on the boundary appear together with $`N`$ minima in the bulk (see Fig. 9a). The minima move towards the center along the line defined by $`\theta =2\pi (k+1/2)/N`$, $`k=1,\mathrm{},N`$, while two saddle points constituting the $`k`$-th pair move symmetrically to the points $`\theta =2\pi k/N`$ and $`\theta =2\pi (k+1)/N`$ as the flux increases in the interval $`\varphi _m(N)<\varphi _e<\varphi _M(N)`$. When $`\varphi _e=\varphi _M(N)`$, with

$$\varphi _M(N)=N\frac{1+x^N}{1x^N}$$
(53)
pairs of saddle points meet at $`\theta =2\pi k/N`$ and for larger values of the flux there are $`N`$ saddle points situated between a vortex and the boundary. In this regime there are $`N+1`$ minima and $`2N`$ saddle points inside the disk as shown in Fig. 9b.

At $`\varphi _e=\varphi _c(N)`$, where $`\varphi _c(N)`$ is given by the formulae (76) and (77), each of the $`N`$ minima coalesce with each of the $`N`$ saddle points at $`r=\rho _c^{1/N},\theta =2\pi (k+1/2)/N`$, where $`\rho _c`$ is given by (76). For $`\varphi _e>\varphi _c(N)`$ there remains only one minimum in the center and $`N`$ saddle points. It can be shown that $`\varphi _c(N)`$ is always greater than $`\varphi _M(N)`$ for $`N>2`$.

Then, for $`\varphi _e=\varphi _c(N)`$ a topological transition occurs, when a minimum and a saddle point coalesce. This transition and the corresponding curve $`\mathrm{\Gamma }`$ is depicted in Fig. 9b and 9c. It is seen that the curve $`\mathrm{\Gamma }`$ undergoes itself a topological transition. For $`\varphi _e<\varphi _c(N)`$ it does not touch the center and has a topology of a ring, while for $`\varphi _e>\varphi _c(N)`$ it has a form of petals of a flower emanating from the center of the disk.

### 5.3 A dynamical interpretation of the curve $`\mathrm{\Gamma }`$

The curve $`\mathrm{\Gamma }`$, defined above as being everywhere normal to the current flow, has a simple interpretation in terms of an effective dynamical system, which gives an alternative method to compute it. Instead of the hamiltonian flow (40), consider the lines of $`\stackrel{}{}h`$, obtained from those of $`\stackrel{}{ȷ}`$ by a rotation of $`90^0`$. This flow is clearly non-Hamiltonian and possess fixed points which are not present in the hamiltonian dynamics, e.g. the positions of the vortices appear in this picture as a point-like sources, though a saddle remain a saddle in this transformation. Then, the curve $`\mathrm{\Gamma }`$ represents a limit cycle of the flow $`\stackrel{}{}h`$ emanating from the saddle point. This is illustrated in Fig. 10 where it is clearly seen that the same curve as in Fig. 7 a) appears as a limit of the lines of $`\stackrel{}{}h`$.

## 6 Conclusion

In this paper we have presented two different aspects of vortex patterns in a mesoscopic superconductor. First we focused on equilibrium properties of a type II superconductor in the London regime. In this case, the Ginzburg-Landau equations reduce to a single linear elliptic equation that can explicitly be solved in simple geometries. We have derived closed expressions for the free energy and the magnetization of a mesoscopic system which have been obtained as limiting cases of a general calculation valid for a system of any size  and not from specific assumptions related to the mesoscopic regime . With these formulae, we have been able to study precisely the matching fields and to compute the paramagnetic Meissner effect, resulting from the existence of metastable states. The second part of the paper has been devoted to the general and topological study of the metastable states using a parallel with the theory of two dimensional dynamical systems. We showed that a vortex pattern can be characterized by as set of topological integers, the number of vortices being one of them. Any change in this set corresponds to a topological transition of the vortex pattern. We calculated the critical fields associated with these transitions and we gave a dual interpretation of these transitions. A key concept is the introduction of a special curve $`\mathrm{\Gamma }`$, that embodies the main geometric features of a configuration of vortices. This curve, which appears mathematically as a limit cycle of the system of currents generated by the vortex pattern, has a simple physical meaning: it separates paramagnetic from diamagnetic domains.

We emphasize that this approach that focuses on topological aspects are generic and can be extended, at least in a qualitative manner, to systems beyond the London regime and with different shapes. Here, we considered a circular cylinder to obtain simple analytical expressions. In recent studies, analytic expressions for the free energy of a cylinder with an elliptic cross-section were found in the mesoscopic limit  and, more generally, an electrostatic analogy was developed that enables to study mesoscopic superconductors in the London regime with the help of conformal transformations ; using such tools, it would be interesting to carry on a quantitative study of topological phase transitions in arbitrary domains.

In this work, only static configurations of vortices have been investigated. A natural extension of the concepts introduced here would be to consider the dynamics of a vortex pattern, where a vortex would move in the field generated by the other vortices. In this context, the study of the deformations of the curve $`\mathrm{\Gamma }`$ would certainly shed some light on the mechanism of vortex nucleation in a superconducting system.

## 7 Acknowledgments

E.A and D.M.G. would like to thank Mark Mineev for fruitful and stimulating discussions. D.M.G. would like to thank SPhT at Saclay for hospitality. This research was supported in part by the U.S.–Israel Binational Science Foundation (BSF), by the Minerva Center for Non-linear Physics of Complex Systems, by the Israel Science Foundation, by the Niedersachsen Ministry of Science (Germany) and by the Fund for Promotion of Research at the Technion.

## Appendix A Magnetic induction in a cylinder

We present the solution of the London equation (15) for the magnetic induction

$$^2h(\stackrel{}{r})2h(\stackrel{}{r})=4\pi \underset{k=1}{\overset{N}{}}\delta (\stackrel{}{r}\stackrel{}{r}_k)$$
(54)
with the boundary condition $`h(R,\theta )=h_e`$. The polar coordinates are convenient for this geometry. The equation (54) with a vanishing right hand side describes the Meissner regime in the absence of vortices. The corresponding solution of the homogeneous equation is given by the $`\theta `$-independent function

$$h_M(r)=h_e\frac{I_0(\sqrt{2}r)}{I_0(\sqrt{2}R)}.$$
(55)
where the functions $`I_n(x)`$, $`K_n(x)`$ are modified Bessel functions of first and second kind of order $`n`$ which provide the regular and singular solutions respectively of (54). The solution of the inhomogeneous equation (54) for $`N`$ vortices is written as a sum of the solution of the homogeneous equation and a particular solution:

$$h(r,\theta )=h_M(r)+h_V(r,\theta )+h_{\overline{V}}(r,\theta )$$
(56)
where $`h_V(r,\theta )`$ is the solution for $`N`$ vortices in the unbound domain

$$h_V(r,\theta )=2\underset{k=1}{\overset{N}{}}K_0(\sqrt{2}|\stackrel{}{r}\stackrel{}{r}_k|)$$
(57)
and the term $`h_{\overline{V}}(r,\theta )`$ is introduced to take care of the boundary conditions:

$$h_{\overline{V}}(R,\theta )=h_V(R,\theta )=2\underset{k=1}{\overset{N}{}}K_0\left(\sqrt{2\left(R^2+r_k^22Rr_k\mathrm{cos}(\theta \theta _k)\right)}\right)$$
(58)
On the other hand, the field $`h_{\overline{V}}`$ can be written as a superposition of solutions of (54) regular in the interior of the boundary:

$$h_{\overline{V}}(R,\theta )=\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}c_nI_n(\sqrt{2}r)e^{in\theta }$$
(59)
where, since the field is real, the relation $`c_n=c_n^{}`$ holds. Expanding the lhs of (58) into Fourier series using standard identities of Bessel functions  we get

$$K_0\left(\sqrt{2\left[R^2+r_k^22Rr_k\mathrm{cos}(\theta \theta _k)\right]}\right)=\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}K_n(\sqrt{2}R)I_n(\sqrt{2}r_k)e^{in(\theta \theta _k)}.$$
(60)
Equating the coefficients of the expansion of both sides of equation (58) one obtains

$$c_n=\frac{K_n(\sqrt{2}R)I_n(\sqrt{2}r_k)}{I_n(\sqrt{2}R)}e^{in\theta _k}$$
(61)
from which the expression (20) follows immediately.

## Appendix B Expression of the Gibbs energy for a cylinder

One takes advantage of the vector identity $`(\times \stackrel{}{h})(\times \stackrel{}{h})=(\stackrel{}{h}\times \times \stackrel{}{h})+\stackrel{}{h}(\times \times \stackrel{}{h})`$ to rewrite $`(\times \stackrel{}{h})(\times \stackrel{}{h})=(\stackrel{}{h}\times \times \stackrel{}{h})\stackrel{}{h}^2\stackrel{}{h}`$, since $`\stackrel{}{h}=0`$. Then one uses the London equation (15) to rewrite the expression (12) for the free energy as

$$=\pi \underset{k=1}{\overset{N}{}}_\mathrm{\Omega }d^2rh(\stackrel{}{r})\delta (\stackrel{}{r}\stackrel{}{r}_k)+\frac{1}{4}_\mathrm{\Omega }d^2r(\stackrel{}{h}\times \times \stackrel{}{h})$$
(62)
The singular part $`_{\mathrm{}}`$ is obtained from the contribution of $`h_V`$ to the first integral. The regular contribution to the free energy can be written using Stokes theorem:

$$_{\mathrm{}}=\pi \underset{k=1}{\overset{N}{}}h_M(\stackrel{}{r}_k)+\pi \underset{k=1}{\overset{N}{}}h_{\overline{V}}(\stackrel{}{r}_k)+\frac{R}{4}_0^{2\pi }𝑑\theta h(R,\theta )h_r(R,\theta )$$
(63)
where $`h_r`$ denotes the partial derivative of $`h(r,\theta )`$ with respect to $`r`$ and the fields $`h_M`$ and $`h_{\overline{V}}`$ are given by (18) and (20). In the last term $`h(R,\theta )=h_e`$ can be taken out of the integration sign and the integral becomes

$`2\pi {\displaystyle \frac{dh_M(R)}{dR}}+2\pi {\displaystyle \frac{}{R}}{\displaystyle _0^{2\pi }}𝑑\theta \left(h_V(R,\theta )+h_{\overline{V}}(R,\theta )\right)`$ (64)
$`=`$ $`2\sqrt{2}\pi h_e{\displaystyle \frac{I_1(\sqrt{2}R)}{I_0(\sqrt{2}R)}}4\sqrt{2}\pi {\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle \frac{I_0(\sqrt{2}r_k)}{I_0(\sqrt{2}R)}}\left(I_0(\sqrt{2}R)K_1(\sqrt{2}R)+K_0(\sqrt{2}R)I_1(\sqrt{2}R)\right)`$
To prove the last equality the properties of the Bessel functions $`I_0^{}(x)=I_1(x)`$, $`K_0^{}(x)=K_1(x)`$ and the expansion (60) were used. Only $`n=0`$ term contributes to the angular integral. Using another property of Bessel functions $`I_0(x)K_1(x)+K_0(x)I_1(x)=1/x`$ and one obtains that:

$$\frac{h_eR}{4}_0^{2\pi }𝑑\theta h_r(R,\theta )=\frac{\sqrt{2}\pi Rh_e^2}{2}\frac{I_1(\sqrt{2}R)}{I_0(\sqrt{2}R)}\pi \underset{k=1}{\overset{N}{}}h_e\frac{I_0(\sqrt{2}r_k)}{I_0(\sqrt{2}R)}$$
(65)
We notice that the contribution from the Meissner field $`h_M`$ in (62) exactly cancels with the second term in the last equation. Hence, the regular part of free energy becomes:

$$_{\mathrm{}}=\sqrt{2}\pi Rh_e^2\frac{I_1(\sqrt{2}R)}{I_0(\sqrt{2}R)}2\pi \underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}\underset{j,k=1}{\overset{N}{}}\frac{K_n(\sqrt{2}R)I_n(\sqrt{2}r_j)I_n(\sqrt{2}r_k)}{I_n(\sqrt{2}R)}\mathrm{cos}n(\theta _j\theta _k)$$
(66)
The term containing the average magnetic induction over the cross section of the cylinder is evaluated along the same lines, using the London equation (15):

$$h_e_\mathrm{\Omega }h(\stackrel{}{r})d^2r=2\pi Nh_e+\frac{h_e}{2}_\mathrm{\Omega }h(\stackrel{}{r})d^2r=\sqrt{2}\pi Rh_e^2\frac{I_1(\sqrt{2}R)}{I_0(\sqrt{2}R)}+2\pi h_e\underset{k=1}{\overset{N}{}}\left(1\frac{I_0(\sqrt{2}r_k)}{I_0(\sqrt{2}R)}\right)$$
(67)
Using the expressions (66) and (67) in the definition (22) of the Gibbs energy completes the derivation of (24).

## Appendix C Study of the critical points

In order to study critical points of the Hamiltonian system defined in section 5.1, one has to solve the equations (45) and (46) for a given number $`N`$ of vortices.

### C.1 Existence and location of the critical points

We first look for critical points which are neither on the boundary of the disk nor at the center, i.e. we concentrate on the generic position $`0<r<1`$. When $`0<r<1`$, one deduces that the equation (46) is satisfied if $`\mathrm{sin}N\theta =0`$, i.e. when $`\theta =k\pi /N`$ for $`k=0,1,\mathrm{},2N1`$. These values of the angle can be subdivided into two groups according to the sign of $`\mathrm{cos}N\theta `$. We now discuss the solution of the equation (45) for the two cases separately.

#### C.1.1 Bulk critical points such that $`\mathrm{cos}N\theta =+1`$.

In this case, we are looking for ’parallel’ critical points that lie in the same direction as the vortices. And we shall prove that for any $`N`$, there exists a critical field $`\varphi _M(N)`$, defined as follows,

$$\varphi _M(N)=N(1+x^N)/(1x^N)$$
(68)
such that when $`\varphi _e>\varphi _M(N)`$ there is always a critical point between each vortex and the boundary, and no such point exists if $`\varphi _e<\varphi _M(N)`$.
Using the new variables $`\xi =x^N`$ and $`\rho =r^N`$ the equation (45) reduces to

$$\frac{\varphi _e}{N}\rho ^{\frac{2}{N}1}=\frac{\xi ^21}{(1\xi \rho )(\xi \rho )}$$
(69)
Looking at this equation, it is clear that systems with $`N=1`$, $`N=2`$ or $`N>2`$ should be discussed separately.

* For $`N=1`$, the plots of the left and right hand sides of this equation are represented in Fig. 11. They cross for $`\rho <1`$ only if the rhs evaluated at $`r=1`$ is less than $`\varphi _e`$, i.e for $`\varphi _e>\varphi _M=(1+\xi )/(1\xi )=(1+x)/(1x)`$. In this case the solution of (69) can be shown to exist using the continuity of the functions. The position of the critical point corresponding to this solution satisfies $`\xi <r<1`$, i.e the point resides between the vortex and the boundary. For $`\varphi _e<\varphi _M`$ there is no solution.
* For $`N2`$, it is worthwhile to consider the reciprocal of (69) which rewrites as

$$\frac{N}{\varphi _e}\rho ^{1\frac{2}{N}}=\frac{1}{\xi ^21}(1\xi \rho )(\xi \rho )$$
(70)
The right hand side is a convex parabola with zeros at $`\xi <1`$ and $`1/\xi >1`$. For $`N=2`$ the lhs is a constant function and it is clear that if $`\varphi _e>\varphi _M=2(1+\xi )/(1\xi )=2(1+x^2)/(1x^2)`$ there is always a solution for $`\xi <\rho <1`$, therefore there exists a critical point between the vortex and the boundary, while $`\varphi _e<\varphi _M`$ no solution exists. For $`N>2`$ the lhs is an increasing fractional power of $`\rho `$ (see Fig. 12). The equation (70) always has a solution between $`\xi `$ and $`1`$ if $`\varphi _e>\varphi _M=N(1+\xi )/(1\xi )=N(1+x^N)/(1x^N)`$ and no solution if $`\varphi _e<\varphi _M`$.

#### C.1.2 Bulk critical points such that $`\mathrm{cos}N\theta =1`$.

We are now looking for ‘anti-parallel’ critical points that lie in a direction bisecting the angle between two neighboring vortices. Here the discussion is more involved:

* when $`N=1`$ there appears a critical field $`\varphi _m(1)=(1x)/(1+x)`$ such that: when $`\varphi _e>\varphi _m(1)`$ there exists a critical point in the direction opposite to the vortex whereas when $`\varphi _e<\varphi _m(1)`$ no such point exists.
* when $`N=2,`$ there are two critical fields $`\varphi _m(2)`$ and $`\varphi _c(2)`$ such that: when $`\varphi _m(2)<\varphi _e<\varphi _c(2)`$ there are exactly two anti-parallel critical points and there are no such points when $`\varphi _e<\varphi _m(2)`$ or $`\varphi _e>\varphi _c(2)`$.
* when $`N>2,`$ one has again two critical fields $`\varphi _m(N)`$ and $`\varphi _c(N)`$ such that: when $`\varphi _e<\varphi _m(N)`$ there exist $`N`$ anti-parallel critical points, when $`\varphi _m(N)<\varphi _e<\varphi _c(N)`$, there are $`2N`$ such points, and when $`\varphi _e>\varphi _c(N)`$, there are no anti-parallel critical points.

We now prove these assertions and calculate the corresponding critical fields.

Consider first that $`N=1`$; then the equation (45) in terms of the variables $`\rho `$ and $`\xi `$ becomes

$$\varphi _e\rho =\frac{1\xi ^2}{(1+\xi \rho )(\xi +\rho )}$$
(71)
And a solution exists in the interval $`0<\rho <1`$ if the value of the lhs at $`\rho =1`$ exceeds the value of the rhs at the same point, as shown in Fig. 13. This happens for $`\varphi _e>\varphi _m(1)=(1x)/(1+x)`$;

Now, when $`N2`$, the equation (45) can be rewritten as

$$\frac{N}{\varphi _e}\rho ^{1\frac{2}{N}}=\frac{1}{1\xi ^2}(1+\xi \rho )(\xi +\rho )$$
(72)
* $`N=2`$: this case is illustrated in Fig. 14. For $`\varphi _e<\varphi _m(2)=2(1\xi )/(1+\xi )=2(1x^2)/(1+x^2)`$ there is no solution, for $`\varphi _e=\varphi _m(2)`$ one solution appears at $`\rho =1`$ and it persists for the range of the flux $`\varphi _m(2)<\varphi _e<\varphi _c(2)=2(1\xi ^2)/\xi =2(1x^4)/x^4`$. For higher flux, $`\varphi _e>\varphi _c(2)`$, this solution disappears at $`\rho =0`$.
* The case $`N>2`$ is the most interesting: Fig. 15 shows the behavior of both sides of the equation (72). Comparing the value of the lhs at $`\rho =1`$ with that of the rhs and using the fact that $`\rho ^{1\frac{2}{N}}`$ has a negative curvature, while that of the rhs is always positive, we can state that for $`\varphi _e<\varphi _m(N)`$, where

$$\varphi _m(N)=N(1\xi )/(1+\xi )=N(1x^N)/(1+x^N)$$
(73)
there is only one solution. Another solution appears at $`\rho =1`$ for $`\varphi _e=\varphi _m(N)`$ and these two solution move towards each other when the flux is increased in the interval $`\varphi _m(N)<\varphi _e<\varphi _c(N)`$, whereas at the critical value of the flux $`\varphi _e=\varphi _c(N)`$ these two solutions coalesce. For $`\varphi _e>\varphi _c(N)`$, the equation (72) has no solutions in the interval $`0<\rho <1`$.

Calculation of the critical flux $`\varphi _c`$:

For $`\varphi _e(N)=\varphi _c`$ the function $`f_{lhs}(\rho )=\frac{N}{\varphi _e}\rho ^{1\frac{2}{N}}`$ and $`f_{rhs}(\rho )=(1+\xi \rho )(\xi +\rho )/(1\xi ^2)`$ are tangent at some point $`\rho =\rho _c`$ as shown in Fig. 15. Equating the value of these functions and the value of their derivatives at $`\rho =\rho _c`$ we obtain the following equations:

$$\{\begin{array}{ccc}\varphi _c(\rho _c+\xi )(1+\xi \rho _c)& =& N(1\xi ^2)\rho _c^{1\frac{2}{N}}\\ \varphi _c(2\xi \rho _c+\xi ^2+1)& =& (N2)(1\xi ^2)\rho _c^{\frac{2}{N}}\end{array}.$$
(74)
Dividing the first equation by the second leads to the quadratic equation for $`\rho _c`$:

$$(N+2)\xi \rho _c^2+2(\xi ^2+1)\rho _c(N2)\xi =0$$
(75)
which has the unique positive solution

$$\rho _c=\frac{(1+\xi ^2)+\sqrt{(1\xi ^2)^2+N^2\xi ^2}}{(N+2)\xi }$$
(76)
The critical flux is obtained by substituting the value of $`\rho _c`$ into the first equation of (74):

$$\varphi _c=\frac{N(1\xi ^2)\rho _c^{1\frac{2}{N}}}{(\rho _c+\xi )(1+\xi \rho _c)}$$
(77)
We notice that for $`N=2`$ the formula (76) yields the correct value $`\rho _c=0`$. Although topologically the cases $`N=2`$ and $`N>2`$ are different, the critical flux $`\varphi _c(2)=(1\xi ^2)/\xi `$, as we shall see later, comes out correctly.

#### C.1.3 Central and boundary critical points

So far we did not discuss the two cases, $`r=0`$ and $`r=1`$. For $`N=1`$ the center, $`r=0`$, of the disk does not satisfy the equations (45),(46) and therefore is not a critical point. For $`N>1`$ the value $`r=0`$ is always a solution of the equations (45),(46). For the points on the boundary, $`r=1`$, the equation (46) is always satisfied and substituting the value $`r=1`$ into (45) we find that there are $`2N`$ critical points on the boundary with an angle given by the following relation:

$$\mathrm{cos}N\theta =\frac{1}{2x^N}\left[1+x^{2N}\frac{N(1x^{2N})}{\varphi _e}\right]$$
(78)
which is meaningful only for values of the flux between $`\varphi _m(N)`$ and $`\varphi _M(N)`$.

### C.2 The nature of the critical points

We now discuss the topological nature of the critical points and classify them by linearizing the stability exponents near each of these points. The problem is then reduced to the study of the Hessian matrix of the second derivatives

$$H=\left(\begin{array}{cc}_r^2h& _r_\theta h\\ _\theta _rh& _\theta ^2h\end{array}\right)$$
(79)
Since at the critical points inside the disk we have $`\mathrm{sin}N\theta =0`$ the mixed derivative vanish and the Hessian matrix is diagonal in polar coordinates. The second derivative with respect to $`\theta `$ is given at the critical points by

$$_\theta ^2h=2N^2\rho \frac{(1\rho ^2)(1\xi ^2)}{(1\xi \rho )^2(\xi \rho )^2}<0$$
(80)
for $`\mathrm{cos}N\theta =+1`$ and by

$$_\theta ^2h=+2N^2\rho \frac{(1\rho ^2)(1\xi ^2)}{(1+\xi \rho )^2(\xi +\rho )^2}>0$$
(81)
for $`\mathrm{cos}N\theta =1`$. The value of the second derivative with respect to $`r`$ can be obtained by looking at the way the first derivative $`_rh`$ changes sign near each solution of the equation (45). We shall now discuss separately the nature of the critical points in the ’parallel’, ’antiparallel’, ’central’ and ’boundary’ cases.

#### C.2.1 Bulk critical points with $`\mathrm{cos}N\theta =+1`$.

For $`\varphi _e>\varphi _M`$, there is only one critical point at $`\rho _1`$ in the interval $`\xi <r<1`$. The results for the second derivatives at this point are summarized in the following table:

$`0<\rho <\rho _1`$ $`\rho _1<\rho <1`$  $`_rh`$  +  $`_r^2h`$ +  $`_\theta ^2h`$   Crit. point Saddle

#### C.2.2 Bulk critical points with $`\mathrm{cos}N\theta =1`$.

The following tables provide a classification of critical points in a variety of regimes.

No critical points.

One critical point at $`\rho _1`$

$`0<\rho <\rho _1`$ $`\rho _1<\rho <1`$  $`_rh`$  +  $`_r^2h`$ +  $`_\theta ^2h`$ +  Crit. point Minimum

No critical points.

One critical point at $`\rho _1`$

$`0<\rho <\rho _1`$ $`\rho _1<\rho <1`$  $`_rh`$  +  $`_r^2h`$ +  $`_\theta ^2h`$ +  Crit. point Minimum

No critical points.

One critical point at $`\rho _1`$

$`0<\rho <\rho _1`$ $`\rho _1<\rho <1`$  $`_rh`$ +   $`_r^2h`$   $`_\theta ^2h`$ +  Crit. point Saddle

Two critical points at $`\rho _1<\rho _2`$

$`0<\rho <\rho _1`$ $`\rho _1<\rho <\rho _2`$ $`\rho _2<\rho <1`$  $`_rh`$ +  +  $`_r^2h`$  +  $`_\theta ^2h`$ + +  Crit. point  Saddle Minimum

No critical points.

#### C.2.3 Central and boundary critical points.

Consider now the point $`r=0`$. Expanding the field $`h(r,\theta )`$ given by (44) for small $`r`$ we obtain

$$h(r,\theta )(\varphi _e2N\mathrm{ln}x)+\varphi _er^2+2\mathrm{cos}(N\theta )\left(\frac{1\xi ^2}{\xi }\right)r^N$$
(82)
As we mentioned before, for $`N=1`$ the point $`r=0`$ is not a critical point (due to the presence of a term linear in $`r`$ ). For $`N=2`$, two last terms are quadratic in $`r`$ and the character of the critical point depends on the value of the flux $`\varphi _e`$. For $`\varphi _e<\varphi _c=(1\xi ^2)/\xi `$ the center, $`r=0`$, is a saddle, while for $`\varphi _e>\varphi _c`$ it turns to be to a minimum. This is to be compared with the prediction of formula 77. It is interesting to notice that this happens when the minimum located along the line $`\mathrm{cos}2\theta =1`$ disappears at the center. For $`N>2`$, the last term in the equation (82) is negligible and the center is always a minimum.

When $`\varphi _m<\varphi _e<\varphi _M`$ we must consider the critical points on the boundary. Since for $`r=1`$ the condition $`_\theta h=0`$ is satisfied for all $`\theta `$, we obtain $`_\theta ^2h=0`$ and $`_r_\theta h0`$ if $`\mathrm{cos}^2N\theta 1`$, which is always satisfied for $`\varphi _e`$ in the interval between $`\varphi _m`$ and $`\varphi _M`$. Thus the determinant of the Hessian matrix (79), given by $`\text{det}H=(_r_\theta h)^2`$, is negative and the eigenvalues are of opposite sign. We conclude that the critical points at the boundary are saddles.


format: text

text:
# CP violating supersymmetric contributions to the electroweak 𝜌 parameter

## Abstract

Effects of CP violation on the supersymmetric electroweak correction to the $`\rho `$ parameter are investigated. To avoid the EDM constraints, we require that arg$`(\mu )<10^2`$ and the non-universal trilinear couplings $`A_f=(0,0,A_0)`$ and also assume that gluinos are heavier than 400 GeV. The CP phase $`\varphi _t=`$ arg($`A_0`$) leads to large enhancement of the relative mass splittings between $`\stackrel{~}{t}_2`$ and $`\stackrel{~}{b}_L(\stackrel{~}{t}_1)`$, which in turn reduces the one-loop contribution of the stop and sbottom to $`\mathrm{\Delta }\rho `$. For small $`\mathrm{tan}\beta `$, such a CP violating effect is prominent. We also study how much the two-loop gluon and gluino contributions are affected by the CP phase. Possible contributions to the $`\rho `$ parameter arising from the Higgs sector with CP violation are discussed.

preprint:     hep-ph/0008073  KIAS-P00036

The minimal supersymmetric standard model (MSSM) is the best motivated extensions of the standard model (SM). As is well known, in the MSSM, there are many new sources of CP violation beyond the Cabibbo-Kobayashi-Maskawa (CKM) phase, arizing from the complex soft supersymmetry (SUSY) breaking terms, i.e., the Majorana gaugino masses $`M_i`$, the trilinear scalar $`A`$ terms, the bilinear scalar $`B`$ terms, as well as from the $`\mu `$ parameter which is the bilinear mixing of the two Higgs chiral superfields in the superpotential. While those supersymmetric CP violating (CPV) phases could give large contributions to the neutron/electron electric dipole moments (EDMs), they are strongly constrained by the current experimental measurements of the neutron/electron EDMs, which for the neutron is $`d_n<1.1\times 10^{25}e`$ cm  and for the electron is $`d_e<4.3\times 10^{27}e`$ cm . To resolve this problem, several possibilities have been suggested. The first one is to make all phases small ($`O(10^2`$)) . But, such small phases would require a significant amount of fine-tuning and are thus undesirable. The second is to take the SUSY spectrum heavy in the several TeV range which lies outside the reach of the accelerators . Another option suggested recently is to use the cancellations among the various contributions to the neutron/electron EDMs, allowing for large CPV phases and a supersymmetric spectrum that is still within the reach of the accelerators . Finally, an interesting possibility is to take a slightly non-universal scenario for the trilinear couplings $`A_f`$ . As shown in Ref. , requiring that arg($`\mu )<10^2`$ and $`A_f=(0,0,A_0)`$, and assuming that gluinos are heavier than about 400 GeV , one may significantly reduce the size of the neutron/electron EDMs due to Weinberg’s three-gluon operator  well below the current experimental limit.

Another CP violation associated with the Higgs boson sector in the MSSM can come from the one-loop corrections of the MSSM Higgs potential through soft CPV Yukawa interactions involving squarks, i.e., $`A`$ terms. As shown in Ref., an immediate consequence of Higgs sector CP violation in the MSSM is the generation of mixing mass terms between the CP even and CP odd Higgs fields. Thanks to those scalar-pseudoscalar mixing mass terms, the small tree-level mass difference between the heaviest Higgs boson and the CP-odd scalar may be lifted considerably. The phenomenological consequences of Higgs sector CP violation have been studied .

A possible way to probe SUSY is to search for the virtual effects of SUSY particles which would be sizable enough to be detected in the present experiments. In most experimentally accessible processes, the heavy SUSY particles decouple from the low energy electroweak observables. However, similar to the SM, if there happens large ”custodial” $`SU(2)_V`$ breaking, the electroweak observables which are parameterized by the $`\rho `$ parameter , defined as the ratio of the strengths of neutral and charged currents at vanishing momentum transfer, can severely be affected. In the MSSM, the main potential source of $`SU(2)_V`$ breaking is the splitting between the stop and sbottom masses. In addition, there can be large splitting within a supersymmetric multiplet, which also affects the $`\rho `$ parameter. From the numerical analysis, it is well known that the leading contribution of the squark loops, in particular stop and sbottom loops, to the $`\rho `$ parameter can reach at the level of a few $`10^3`$ which lies within the range of the experimental sensitivity . There are also small additional contributions coming from the exchange of the additional Higgs bosons, the charginos and neutralinos. Some discussions of supersymmetric contributions to the $`\rho `$ parameter already exist in the literature . Moreover, higher order corrections have also been studied so as to treat the SUSY loop contributions to the electroweak observables at the same level of accuracy as the SM contribution . However, so far, the analyses have been done within the context of the CP conserving (CPC) MSSM. Since new CPV phases may affect the SUSY spectrum and generate additional couplings which do not appear in the CPC MSSM, it would be quite interesting to study the effects of the CPV SUSY phases on the electroweak observables via the $`\rho `$ parameter.

In this respect, the purpose of this letter is to examine how much the $`\rho `$ parameter can be affected by the non-vanishing CP phases in the MSSM. In this analysis, we will impose the universal conditions on the soft SUSY breaking terms which provide three complex parameters, i.e., universal gaugino mass $`m_{1/2}`$, the trilinear soft term $`A_0`$ and the bilinear soft term $`B_0`$. In addition, there is a complex mass parameter $`\mu `$. However, not all phases of the four complex parameters are physical . It is always possible to make a phase transformation on the gaugino fields so as to make $`m_{1/2}`$ real. From minimization conditions on the Higgs potential, one can make the phase of $`B_0`$ equal to the phase difference of the two Higgs doublets in the MSSM. As a result, there are only two independent CPV phases in the MSSM with universal soft SUSY breaking terms, which are usually chosen to be arg($`A_0`$) and arg($`\mu `$). To satisfy the EDM constraints, we will simply take arg($`\mu `$) to be zero and the trilinear coupling of the Higgs bosons to the squarks of the first and second generation to be much smaller than the one of the third generation, as mentioned above, and assume that the gluino mass to be larger than 400 GeV. Thanks to this non-vanishing phase arg($`A_0`$), as will be shown later, the stop spectrum and their mixing angle are changed, which lead in turn to the shift of the contributions to $`\mathrm{\Delta }\rho `$. Also, since the neutral Higgs boson mass eigenstates in the CPV MSSM are mixtures of CP even and CP odd states, new gauge-Higgs couplings are induced at the tree level through those mixings. These couplings can generate the additional radiative corrections to the gauge boson self-energies, and may thus affect the $`\rho `$ parameter. However, the contributions coming from the gaugino and higgsino exchanges will not be affected by the non-vanishing CPV phase arg($`A_0`$) in the leading order because they are affected by only the CPV phase arg($`\mu `$) which is neglected in this analysis. Thus, we shall not consider those contributions. In addition, we shall study how much the dominant two-loop contributions via gluon and gluino exchanges can be shifted by the effect of the CPV phases. The numerical calculation will be done with the help of the program package FEYNHIGGS  which is modified so as to calculate the CPV SUSY contributions to the $`\rho `$ parameter.

The $`\rho `$ parameter is represented by

$$\rho =\frac{1}{1\mathrm{\Delta }\rho },\mathrm{\Delta }\rho =\frac{\mathrm{\Pi }_{ZZ}(0)}{M_Z^2}\frac{\mathrm{\Pi }_{WW}(0)}{M_W^2},$$
(1)
where $`\mathrm{\Pi }_{WW(ZZ)}(0)`$ denotes the transverse parts of the $`W/Z`$-boson self-energies at zero momentum transfer. In the SM, the main contribution comes from the large splitting between the top quark and the bottom quark masses and is given by

$$\mathrm{\Delta }\rho _0^{SM}=\frac{N_cG_F}{8\sqrt{2}\pi ^2}F_0(m_t^2,m_b^2),$$
(2)
where $`N_c`$ is the color factor and the function $`F_0`$ is given by

$$F_0(x,y)=x+y\frac{2xy}{xy}\mathrm{ln}\frac{x}{y}.$$
(3)
In the limit of large mass splitting it becomes proportional to the heavy quark mass squared, i.e., $`F_0(m_t^2,m_b^2)m_t^2`$. Including QCD corrections, the SM contribution of $`\mathrm{\Delta }\rho `$ is $`\mathrm{\Delta }\rho _1^{SM}=\mathrm{\Delta }\rho _0^{SM}\left(\frac{2\alpha _s}{3\pi }\right)(1+\pi ^2/3)`$ . For $`\alpha _s(M_Z)0.12`$, the QCD correction reduces the one-loop result by about $`10\%`$ and shifts $`m_t`$ upwards by about $`10`$ GeV .

(A) Squark contributions: Since the scalar partners of the light quarks are expected to be almost degenerate in most supersymmetric theories, their contributions to $`\mathrm{\Delta }\rho `$ are negligible and thus only third generation will contribute. The stop mass matrix is given by

$`M_{\stackrel{~}{t}}^2=\left(\begin{array}{cc}m_{\stackrel{~}{t}_L}^2+m_t^2+D_L& m_tm_{LR}^{}\\ m_tm_{LR}& m_{\stackrel{~}{t}_R}^2+m_t^2+D_R\end{array}\right)`$ (4)
where $`m_{LR}=A_0\mu ^{}/\mathrm{tan}\beta `$, $`D_L=\left(\frac{1}{2}\frac{2}{3}\mathrm{sin}^2\theta _W\right)\mathrm{cos}2\beta M_Z^2`$ and $`D_R=\frac{2}{3}\mathrm{sin}^2\theta _W\mathrm{cos}2\beta M_Z^2`$. The stop fields $`(\stackrel{~}{t}_L,\stackrel{~}{t}_R)`$ are linear combinations of the mass eigenstates $`\stackrel{~}{t}_{1,2}`$ so that $`\stackrel{~}{t}_L=\mathrm{cos}\theta _{\stackrel{~}{t}}\stackrel{~}{t}_1\mathrm{sin}\theta _{\stackrel{~}{t}}e^{i\phi _t}\stackrel{~}{t}_2,\stackrel{~}{t}_R=\mathrm{sin}\theta _{\stackrel{~}{t}}e^{i\phi _t}\stackrel{~}{t}_1+\mathrm{cos}\theta _{\stackrel{~}{t}}\stackrel{~}{t}_2`$, where

$$\mathrm{tan}2\theta _{\stackrel{~}{t}}=\frac{2m_t|m_{LR}|}{m_{\stackrel{~}{t}_L}^2m_{\stackrel{~}{t}_R}^2+D_LD_R},$$
(5)
and the mass eigenvalues of stop fields are given by

$`m_{\stackrel{~}{t}_{1,2}}^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}[m_{\stackrel{~}{t}_L}^2+m_{\stackrel{~}{t}_R}^2+2m_t^2+D_L+D_R`$ (7)
$`([m_{\stackrel{~}{t}_L}^2m_{\stackrel{~}{t}_R}^2+D_LD_R]^2+4m_t^2|m_{LR}|^2)^{1/2}],`$
$`|m_{LR}|`$ $`=`$ $`|A_0\mu ^{}/\mathrm{tan}\beta |`$ (8)
$`=`$ $`\sqrt{|A_0|^2+|\mu |^2/\mathrm{tan}^2\beta 2|A_0||\mu |/\mathrm{tan}\beta \mathrm{cos}\varphi _t},`$ (9)
where $`\varphi _t\text{arg}(A_0)`$. As expected, the non-vanishing phase $`\varphi _t`$ will change the stop mass spectrum which can in turn lead to the shift of $`\mathrm{\Delta }\rho `$. Neglecting the mixing in $`\stackrel{~}{b}`$ sector, the contribution of $`(\stackrel{~}{t},\stackrel{~}{b})`$ to $`\mathrm{\Delta }\rho `$ is presented at one-loop order by:

$`\mathrm{\Delta }\rho _{\stackrel{~}{t}\stackrel{~}{b}}`$ $`=`$ $`{\displaystyle \frac{3G_F}{8\sqrt{2}\pi ^2}}[\mathrm{sin}^2\theta _{\stackrel{~}{t}}\mathrm{cos}^2\theta _{\stackrel{~}{t}}F_0(m_{\stackrel{~}{t_1}}^2,m_{\stackrel{~}{t_2}}^2)`$ (10)
$`+`$ $`\mathrm{cos}^2\theta _{\stackrel{~}{t}}F_0(m_{\stackrel{~}{t_1}}^2,m_{\stackrel{~}{b_L}}^2)+\mathrm{sin}^2\theta _{\stackrel{~}{t}}F_0(m_{\stackrel{~}{t_2}}^2,m_{\stackrel{~}{b_L}}^2)],`$ (11)
Before examining how much the value of $`\mathrm{\Delta }\rho _{\stackrel{~}{t}\stackrel{~}{b}}`$ can be affected by the phase $`\varphi _t`$, let us first investigate the mass splittings between the two squark mass eigenstates which may indicate to what extent the contribution of $`(\stackrel{~}{t},\stackrel{~}{b})`$ to $`\mathrm{\Delta }\rho `$ is deviated. In Fig. 1 a (b), we present numerical estimation of the relative mass splittings between the two squark mass eigenstates, $`|m_{\stackrel{~}{t}_1}^2m_{\stackrel{~}{t}_2}^2|/m_{\stackrel{~}{q}}^2`$ and $`|m_{\stackrel{~}{t}_2}^2m_{\stackrel{~}{b}_L}^2|/m_{\stackrel{~}{q}}^2`$, as a function of the phase $`\varphi _t`$, for $`\mathrm{tan}\beta =1.6(20),m_{\stackrel{~}{q}}=200`$ GeV, $`|A_0|=200`$ GeV and $`|\mu |=100`$ GeV. Here we take $`0\varphi _t\pi `$. As the phase $`\varphi _t`$ increases, the both relative mass splittings increase. The effect of the maximal CP violation corresponding to $`\varphi _t=\pi `$ may lead to the enhancement of the relative mass splittings by about $`90\%(40\%)`$ for $`|m_{\stackrel{~}{t}_1}^2m_{\stackrel{~}{t}_2}^2|/m_{\stackrel{~}{q}}^2`$ ($`|m_{\stackrel{~}{t}_2}^2m_{\stackrel{~}{b}_L}^2|/m_{\stackrel{~}{q}}^2`$) compared to the CPC case. This is attributed to the fact that the absolute value of $`m_{LR}`$ monotonically increases with the increasing phase $`\varphi _t`$, which makes the mass splitting between $`\stackrel{~}{t}_1(\stackrel{~}{b}_L)`$ and $`\stackrel{~}{t}_2`$ larger, and the mixing angle $`\theta _{\stackrel{~}{t}}`$ closer to $`\pi /4`$. Note that the mixing angle $`\theta _{\stackrel{~}{t}}`$ close to $`\pi /4`$ are naturally obtained because the off-diagonal elements of the stop mass matrix may be larger than the difference of the diagonal entries. As one can see from Fig. 1 a, in the case of small $`\mathrm{tan}\beta `$, the splitting $`|m_{\stackrel{~}{t}_2}^2m_{\stackrel{~}{b}_L}^2|/m_{\stackrel{~}{q}}^2`$ is larger than $`|m_{\stackrel{~}{t}_1}^2m_{\stackrel{~}{t}_2}^2|/m_{\stackrel{~}{q}}^2`$ for $`\varphi _t=0`$ (CPC case), whereas the latter becomes larger than the former for large $`\varphi _t`$. Thus, the splitting originating from the left-right mixing of stops may be more important for the CPV case with large phase. In the case of large $`\mathrm{tan}\beta `$, $`|m_{\stackrel{~}{t}_1}^2m_{\stackrel{~}{t}_2}^2|/m_{\stackrel{~}{q}}^2`$ is larger than $`|m_{\stackrel{~}{t}_2}^2m_{\stackrel{~}{b}_L}^2|/m_{\stackrel{~}{q}}^2`$ for all $`\varphi _t`$, and the mass splittings are weakly dependent on $`\varphi _t`$ since the term concerned with the phase $`\varphi _t`$ in $`|m_{LR}|`$ is suppressed by the large $`\mathrm{tan}\beta `$. We also observed that the mass splitting $`|m_{\stackrel{~}{t}_1}^2m_{\stackrel{~}{b}_L}|/m_{\stackrel{~}{q}}^2`$ is much smaller than the other two. For the same input values of the SUSY parameters, the dependence of $`\mathrm{\Delta }\rho _{\stackrel{~}{t}\stackrel{~}{b}}`$ on the phase $`\varphi _t`$ is shown as solid lines in Fig. 2. Although the CPV case provides larger mass splitting than the CPC case, $`\mathrm{\Delta }\rho _{\stackrel{~}{t}\stackrel{~}{b}}`$ is monotonically decreased with the increasing $`\varphi _t`$. The effect of the CPV phase reduces the value of $`\mathrm{\Delta }\rho _{\stackrel{~}{t}\stackrel{~}{b}}`$ for the CPC case by about up to $`35\%`$. The reason is that as the phase increases, the mixing angle $`\theta _{\stackrel{~}{t}}`$ gets close to the maximal mixing and $`F_0(m_{\stackrel{~}{t}_1}^2,m_{\stackrel{~}{t}_2}^2)`$ increases faster than $`F_0(m_{\stackrel{~}{t}_2}^2,m_{\stackrel{~}{b}_L}^2)`$ in Eq. (8) so that they lead to destructive contribution to $`\mathrm{\Delta }\rho _{\stackrel{~}{t}\stackrel{~}{b}}`$. As $`|A_0|`$ increases, $`\mathrm{\Delta }\rho _{\stackrel{~}{t}\stackrel{~}{b}}`$ tends to decrease, whereas it tends to increase with the increasing $`|\mu |`$. We have observed that the effect of the CPV phase is diminished as long as $`m_{\stackrel{~}{q}}`$ becomes larger than $`|A_0|`$ and $`|\mu |`$. This is because $`m_{\stackrel{~}{t}_1}m_{\stackrel{~}{t}_2}m_{\stackrel{~}{b}_L}`$ in the limit of large $`m_{\stackrel{~}{q}}`$.

The two-loop QCD correction to $`\mathrm{\Delta }\rho _{\stackrel{~}{t}\stackrel{~}{b}}`$ is given by Eq.(8) in Ref., and the two-loop contribution mediated by gluino exchange is also represented in Ref.. Those gluon and gluino contributions add up to about $`30\%`$ of the one-loop value in the CPC case.

In Fig. 3, the dependence of the dominant two-loop QCD correction and gluino contribution to the $`\rho `$ parameter on the phase $`\varphi _t`$ is plotted by the dashed line and dotted line, respectively, for $`\mathrm{tan}\beta =1.6`$ (a) and $`\mathrm{tan}\beta =20`$ (b). As the phase $`\varphi _t`$ increases, the two-loop gluonic SUSY contribution decreases, whereas the gluino contribution increases. The reason that the gluino contribution increases with the increasing phase $`\varphi _t`$ is that its contribution depends inversely on the mass splittings between $`\stackrel{~}{t}_1(\stackrel{~}{b}_L)`$ and $`\stackrel{~}{t}_2`$. In particular, for small $`\mathrm{tan}\beta `$, the value of the gluino contribution becomes larger than that of the gluonic SUSY contribution at large $`\varphi _t`$.

(B) Higgs sector contributions: There are also possible contributions arising from the Higgs sector to the $`\rho `$ parameter. In the CPC case, it is known  that the Higgs boson masses and couplings to gauge bosons are related in such a way that they lead to large cancellations in their contributions to the $`\rho `$ parameter, which is at most of order of $`10^4`$. In the decoupling limit where the heavy neutral CP-even Higgs, CP-odd Higgs and the charged Higgs bosons are nearly degenerate and their couplings to gauge bosons tend to zero, while the lightest CP-even Higgs boson reaches its maximal mass value and has almost the same properties as the SM Higgs boson. Then, the contribution of the Higgs sector of the MSSM to the $`\rho `$ parameter is practically the same as in the SM. However, while the contribution of the SM Higgs sector gives rise to logarithmic $`\mathrm{log}M_h/M_Z`$ in the limit of large Higgs mass which may reach about $`10^3`$ order, the MSSM contribution reaches at most $`10^4`$ due to the upper bound on the mass of the lightest Higgs boson.

Now, let us take into account the contributions coming from the neutral Higgs sector with CP violation. The decomposition of the two Higgs doublets is given by

$`\mathrm{\Phi }_1=\left(\begin{array}{c}\varphi _1^+\\ \frac{1}{\sqrt{2}}(v_1+\varphi _1+ia_1)\end{array}\right),\mathrm{\Phi }_2=e^{i\theta }\left(\begin{array}{c}\varphi _2^+\\ \frac{1}{\sqrt{2}}(v_2+\varphi _2+ia_2)\end{array}\right),`$ (12)
where $`v_1`$ and $`v_2`$ are the vacuum expectation values and $`\theta `$ is their relative phase. As shown in Ref. , the Higgs quartic couplings receive significant radiative corrections from $`A_f`$ terms. The vacuum expectation values $`v_1,v_2`$ and the phase $`\theta `$ are determined from the minimization conditions on the Higgs potential of the MSSM. The CP odd fields are rotated to $`a_1=\mathrm{cos}\beta G^0\mathrm{sin}\beta a`$ and $`a_2=\mathrm{sin}\beta G^0+\mathrm{cos}\beta a`$ so that the Higgs potential shows up a flat direction with respect to the Goldstone field $`G^0`$. The non-vanishing CPV phases in $`A_t`$ and $`A_b`$ can lead to the non-vanishing $`\theta `$, for which the scalar-pseudoscalar mixing mass terms are generated and the neutral Higgs boson mass matrix in the weak basis $`(G^0,a,\varphi _1,\varphi _2)`$ can be written by

$`M_0^2=\left(\begin{array}{cc}M_P^2& M_{PS}^2\\ M_{SP}^2& M_S^2\end{array}\right),`$ (13)
where $`M_P^2`$ and $`M_S^2`$ are $`2\times 2`$ matrix for CPC parts in the basis $`(G^0,a)`$ and $`(\varphi _1,\varphi _2)`$, respectively, and $`M_{PS}^2=(M_{SP}^2)^T`$ describes the CPV mixings between $`(G^0,a)`$ and $`(\varphi _1,\varphi _2)`$. From the tadpole condition, the Goldstone field $`G^0`$ does not mix with the other neutral fields and thus becomes massless. Then, the neutral Higgs mass matrix $`M_0^2`$ reduces to a $`3\times 3`$ matrix $`M^2`$, which is spanned in the basis $`(a,\varphi _1,\varphi _2)`$. The Higgs boson mass matrix $`M^2`$ can be diagonalized with the help of an orthogonal rotation matrix $`U`$: $`U^TM^2U=diag(M_{H_3}^2,M_{H_2}^2,M_{H_1}^2)`$.

Then, the contributions of the Higgs sector with CP violation to the $`\rho `$ parameter can be given by

$`\mathrm{\Delta }\rho _H`$ $`=`$ $`{\displaystyle \frac{3\alpha }{16\pi c_W^2}}{\displaystyle \underset{i=1}{\overset{3}{}}}U_{1i}^2\mathrm{\Delta }\rho _H^{SM}(M_{H_i}^2)+{\displaystyle \frac{\alpha }{16\pi s_W^2M_W^2}}`$ (14)
$`\times `$ $`\left({\displaystyle \underset{i=1}{\overset{3}{}}}(1U_{1i}^2)F_0(M_{H_i}^2,M_{H_c}^2){\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{\overset{3}{}}}(U_{1i}U_{2j}U_{2i}U_{1j})^2F_0(M_{H_i}^2,M_{H_j}^2)\right)`$ (15)
where $`M_{H_c}`$ denotes the charged Higgs boson mass and

$`\mathrm{\Delta }\rho _H^{SM}(M^2)`$ $`=`$ $`{\displaystyle \frac{\alpha }{16\pi s_W^2M_W^2}}[F_0(M^2,M_W^2)F_0(M^2,M_Z^2)]`$ (16)
$`+`$ $`{\displaystyle \frac{\alpha }{16\pi s_W^2}}\left[{\displaystyle \frac{M^2}{M^2M_W^2}}\mathrm{log}{\displaystyle \frac{M^2}{M_W^2}}{\displaystyle \frac{1}{c_W^2}}{\displaystyle \frac{M^2}{M^2M_Z^2}}\mathrm{log}{\displaystyle \frac{M^2}{M_Z^2}}\right]`$ (17)
corresponds to the SM Higgs boson contribution to $`\mathrm{\Delta }\rho `$.

In Fig. 3, we present the dependence of the Higgs boson sector contribution to $`\mathrm{\Delta }\rho `$ on the phase $`\varphi _t`$ (solid lines). Similar to the one-loop contribution of stop quark to $`\mathrm{\Delta }\rho `$, the one-loop contribution of Higgs sector is decreased with the increasing phase $`\varphi _t`$. This is mainly because the relative mass splittings among three neutral Higgs bosons decrease as the phase $`\varphi _t`$ increases as shown in Fig. 4. In the case of small $`\mathrm{tan}\beta `$, $`|M_3^2M_2^2|/M_{H_c}^2`$ (solid line) is larger than $`|M_2^2M_1^2|/M_{H_c}^2`$ (dotted line). For only small $`\mathrm{tan}\beta `$, the contribution of the Higgs sector is larger than the two-loop contributions mediated by the gluon and gluino exchanges, and those three contributions are almost the same at $`\varphi _t=\pi `$, which increase the one-loop contribution of stop and sbottom by about $`3035\%`$. For large $`\mathrm{tan}\beta `$, the Higgs sector contribution to $`\mathrm{\Delta }\rho `$ is negligibly small.

In summary, effects of CP violation on the supersymmetric electroweak correction to the $`\rho `$ parameter have been investigated. To avoid the EDM constraints, we have required negligibly small arg$`(\mu )`$ and the non-universal couplings $`A_f=(0,0,A_0)`$ and also assumed that the mass of gluino is larger than 400 GeV. The CP phase $`\varphi _t`$ leads to large enhancement of the relative mass splittings between $`\stackrel{~}{t}_2`$ and $`\stackrel{~}{b}_L(\stackrel{~}{t}_1)`$, which in turn reduces the one-loop contribution of the stop and sbottom to the $`\rho `$ parameter. For small $`\mathrm{tan}\beta `$, such a CP violating effect is prominent. We have also studied how much the two-loop gluon and gluino contributions are affected by the CP phase. Possible contributions to the $`\rho `$ parameter arising from the Higgs sector with CP violation have been discussed.

The work of JDK is supported in part by the Korea Institute for Advanced Study.