text:
# Evolution of highly buoyant thermals in a stratified layer

## 1 Introduction

Highly buoyant bubbles with large specific entropy excess relative to the surroundings have been invoked by Hazlehurst (1985) in an attempt to explain the almost equal effective temperatures of the two components of contact binaries.

As noted by Sinjab et al. (1990) there are strong parallels between the highly buoyant bubbles of Hazlehurst and the ‘thermals’ found to occur in the earth’s atmosphere. Subsequently, Hazlehurst (1990) confirmed the existence of a formal relationship between Turner’s (1963) treatment of thermals, involving entrainment of matter, and his own treatment of highly buoyant bubbles (‘interlopers’) in contact binaries, as these bubbles annex new material.

In this paper we shall discuss the question of the validity, from a fluid-dynamical standpoint, of the simple bubble or thermal picture. We shall then go on to show how it is possible to determine numerically the value of the entrainment coefficient $`\alpha `$ (here called $`\alpha _\mathrm{v}`$) which enters Turner’s and several other investigations; the determination of a related coefficient (here called $`\alpha _\mathrm{m}`$) entering the Hazlehurst theory is also discussed.

We believe this to be the first attempt to evaluate the (previously semi-empirical) entrainment coefficient $`\alpha `$ on a fluid-dynamical basis.

## 2 Model setup

We adopt a basic setup of our model that is similar to that used normally to study convection in a stratified plane-parallel layer between impenetrable boundaries (e.g. Hurlburt et al. 1984). In particular, we use stress-free boundary conditions at the top and bottom, with a prescribed flux $`F`$ at the bottom and a prescribed temperature $`T_{\mathrm{top}}`$ at the top. Here, however, we assume the thermal equilibrium stratification to be adiabatic, so it is marginally stable to the onset of convection. Thus, when we insert a hot (buoyant) bubble, it will rise unaffected by the stratification (except for effects related to the growth and expansion of the ascending bubble), so there is no restoring force acting on the bubble as it rises.

### 2.1 Adiabatic stratification

Hydrostatic equilibrium requires that the weight of the atmosphere is balanced by the pressure gradient. However, if the entropy is constant, the pressure gradient can be written as $`\rho ^1\mathbf{}p=\mathbf{}h`$, where $`p`$ is pressure, $`\rho `$ is density, and $`h`$ is enthalpy. We adopt a perfect gas for which $`h=c_\mathrm{p}T`$, where $`c_\mathrm{p}`$ is the specific heat at constant pressure. For constant gravitational acceleration $`g>0`$, this implies a linear temperature stratification that is given by

$$\frac{\mathrm{d}}{\mathrm{d}z}(c_\mathrm{p}T)=g,$$
(1)
where $`z`$ is depth, which is assumed to increase downwards. Thus, the vertical temperature profile of the basic state is given by

$$T=T_{\mathrm{top}}+(zz_{\mathrm{top}})g/c_\mathrm{p}.$$
(2)
In the absence of any motion there is only radiative flux, $`𝑭`$, for which we adopt the diffusion approximation, so

$$𝑭=K\mathbf{}T,$$
(3)
where $`K`$ is the radiative conductivity. Thermal equilibrium requires $`\mathbf{}𝑭=0`$, so the $`z`$-component of the flux is constant but, because the temperature gradient is constant, this is only possible if $`K=\mathrm{const}`$. We now need to discuss the choice of $`K`$ and other parameters.

### 2.2 Choice of parameters

We adopt nondimensional units by defining a unit length $`d`$, and our bubble will usually have the initial radius $`R_0=d`$ (although we present initially some cases where $`R_0=0.5d`$. We measure time in units of $`(d/g)^{1/2}`$, density in units of $`\rho _0`$ (we choose $`\rho =\rho _0`$ at the location where the centre of the bubble will be introduced), and specific entropy, $`s`$, in units of $`c_\mathrm{p}`$. This corresponds to setting

$$d=g=\rho _0=c_\mathrm{p}=1.$$
(4)
A relevant nondimensional quantity is the normalised input flux

$$\frac{F}{\rho _0(gd)^{1/2}}\text{(nondimensional input flux)}.$$
(5)
For secondaries of contact binaries this ratio is around $`10^3`$. Specifying $`F`$ fixes $`K=c_\mathrm{p}F/g`$. There is however the numerical constraint that the mesh Peclet number, based on the sound speed $`c_\mathrm{s}`$ and the mesh width $`\mathrm{\Delta }x`$,

$$\text{Pe}_{\mathrm{grid}}=\mathrm{\Delta }xc_\mathrm{s}/\chi ,$$
(6)
should not exceed a certain empirical upper limit of 10–100. This limit arises from the fact that the advection of sharp structures tends to generate small ripples on the scale of the mesh which need to be damped. Here, $`c_\mathrm{s}^2=\gamma p/\rho `$, $`\gamma =c_\mathrm{p}/c_\mathrm{v}`$ is the ratio of specific heats, and $`\chi =K/\rho c_\mathrm{p}`$ is the radiative diffusivity. In the present paper we shall consider values of $`F`$ in the range 0.001–0.005. In all cases we assume $`\gamma =5/3`$.

We adopt cartesian coordinates $`(x,y,z)`$ where $`z`$ points downwards. The bubble centre is placed initially at $`x=y=z=0`$. In order that the bubble be initially sufficiently far away from the bottom boundary and that it can rise over a distance of at least a few times its own radius (which is unity in most of the models), we chose the vertical extent of the computational box to be from $`z_{\mathrm{top}}=6`$ to $`z_{\mathrm{bot}}=2`$.

Next, we fix the temperature stratification within the box by specifying the value of $`T_{\mathrm{top}}`$, which is held constant by the boundary condition chosen. The choice of this parameter is restricted by numerical considerations. Since temperature is proportional to the pressure scale height, small values of $`T_{\mathrm{top}}`$ imply a short pressure scale height at the surface. However, numerical stability and accuracy considerations require that the scale height cannot be less than at least 2–3 mesh zones. The non-dimensional pressure scale height at the top is

$$\xi _0=H_{\mathrm{p},\mathrm{top}}/d=(c_\mathrm{p}c_\mathrm{v})T_{\mathrm{top}}/gd,$$
(7)
where we have included $`d`$ and $`g`$ factors, even though they are unity. For runs with $`N_z=50`$ meshpoints in the vertical direction, we were able to use $`\xi _0=0.3`$. For orientation, we note that this yields a $`\mathrm{ln}(p_{\mathrm{bot}}/p_{\mathrm{top}})`$ of 6.2 pressure scale heights between top and bottom of the box. The local pressure scale height at $`z=0`$ is then given by $`0.4\times |z_{\mathrm{top}}|+\xi _0=2.7`$. Using Eqs (2) and (7) we also see that the top of the adiabatic atmosphere would be at

$$z_{\mathrm{}}=z_{\mathrm{top}}2.5\xi _0\text{(for }\gamma =5/3\text{)},$$
(8)
but this should be outside the computational box. For $`z_{\mathrm{top}}=6`$ and $`\xi _0=0.3`$ this gives $`z_{\mathrm{}}=6.75`$.

For a perfect gas we have $`p/\rho =(c_\mathrm{p}c_\mathrm{v})T`$ and $`s=c_\mathrm{v}\mathrm{ln}pc_\mathrm{p}\mathrm{ln}\rho `$ (in dimensional form). The density stratification follows then from hydrostatic balance and the assumption that the entropy of the unperturbed model is constant:

$$\rho /\rho _0=(1z/z_{\mathrm{}})^{1/(\gamma 1)}.$$
(9)
We mentioned already that we express density in units of $`\rho _0`$, which is the density at $`z=0`$, which is where the bubble will be positioned. Using this together with the definition of $`s`$ we can use Eq. (8) to express the value of the background entropy,

$$s_0=0.6\mathrm{ln}(0.4z_{\mathrm{}})\text{(for }\gamma =5/3\text{)}.$$
(10)
For $`z_{\mathrm{}}=6.75`$ this gives $`s_0=0.596`$.

This completes the definition of the background state. We now turn to the discussion of the model equations and initial and boundary conditions adopted.

### 2.3 Governing equations

In the dynamical case the specific entropy is not only affected by the radiative flux divergence, but also by the rate of viscous dissipation, so

$$\rho T\frac{\mathrm{D}s}{\mathrm{D}t}=\mathbf{}K\mathbf{}T+2\nu \rho 𝗦{}_{}{}^{2},$$
(11)
where $`\mathrm{D}/\mathrm{D}t=/t+𝒖\mathbf{}`$ is the lagrangian derivative, $`\nu =\text{const}`$ is the kinematic viscosity and

$$𝖲_{ij}=\frac{1}{2}(_ju_i+_iu_j\frac{2}{3}\delta _{ij}_ku_k)$$
(12)
is the (traceless) rate of strain tensor. Equation (11) is solved together with the momentum equation,

$$\frac{\mathrm{D}𝒖}{\mathrm{D}t}=\frac{1}{\rho }\mathbf{}p+𝒈+\frac{1}{\rho }\mathbf{}(2\nu \rho 𝗦),$$
(13)
and the continuity equation

$$\frac{\mathrm{D}\mathrm{ln}\rho }{\mathrm{D}t}=\mathbf{}𝒖.$$
(14)
We solve Eqs (11), (13) and (14) using the sixth order compact derivative scheme of Lele (1992) and a third order Hyman scheme for the time step. For earlier applications of this code see Nordlund & Stein (1990) and Brandenburg et al. (1996).

The value of $`\nu `$ is dictated again by numerical considerations, and in practice we take $`\nu \chi `$, but note that $`\nu `$ is independent of $`z`$ whilst $`\chi `$ is not. In all cases considered below we have used $`\nu =6\times 10^3`$.

### 2.4 Initial and boundary conditions

The bubble centre is placed initially at $`x=y=z=0`$ and has an entropy profile of the form

$$s=s_0+\mathrm{\Delta }s\mathrm{\Phi }(r),$$
(15)
where

$`\mathrm{\Phi }(r)=\{\begin{array}{cc}\mathrm{exp}\left({\displaystyle \frac{r^2}{R_0^2r^2}}\right)\hfill & \text{if }r<R_0\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}`$ (18)
is a profile function defining the initial shape of the bubble. Here $`r`$ is the initial distance from the centre. Initial pressure equilibrium requires that the increase of $`s`$ is compensated by a corresponding decrease of $`\mathrm{ln}\rho `$, so

$$\mathrm{ln}\rho =\frac{1}{\gamma 1}\mathrm{ln}\left(1\frac{z}{z_{\mathrm{}}}\right)\mathrm{\Delta }s\mathrm{\Phi }(r).$$
(19)
The initial entropy excess of the bubble is given by the parameter $`\mathrm{\Delta }s`$. \[In all cases presented we take $`\mathrm{\Delta }s=0.5`$, which corresponds to the value used by Hazlehurst (1985), who adopts units where the nondimensional specific entropy is larger by a factor of 5. We note that larger values of $`\mathrm{\Delta }s`$ make the bubble rise faster, but we found that even for $`\mathrm{\Delta }s=2`$ the motion remained subsonic.\]

We assume stress-free boundaries at the top and bottom, and prescribe the value of $`T`$ at the top and $`T/z`$ at the bottom, as is usual for convection calculations (e.g. Hurlburt et al. 1984). We use periodic boundary conditions in the $`x`$ and $`y`$ directions. The horizontal extent of the box, $`|x|<L_x`$ and $`|y|<L_y`$, is varied between $`L_x=L_y=4`$ and 16.

### 2.5 Allowing for three-dimensional effects

In order to assess the fragility of the bubble during its ascent we have adopted in many cases substantial initial velocity perturbations. In Fig. 1 we show a three-dimensional representation of the entropy for a run with $`\mathrm{\Delta }s=0.5`$ using initial velocity perturbations with $`u_{\mathrm{max}}/c_\mathrm{s}=0.4`$ and $`u_{\mathrm{rms}}/c_\mathrm{s}=0.02`$. As is evident from Fig. 1, the entropy of the blob is hardly affected by these perturbations and only near the surface does one see strong perturbations.

It is interesting that the bubble remains an entity during much of its ascent. In fact, even when the initial condition is quite different, bubble-like structures tend to develop. As an example we show in Fig. 2 a case where we have introduced an almost uniform horizontal layer of enhanced specific entropy with a gaussian vertical profile initially. We have superimposed random small scale perturbations to get the buoyancy instability started.

In the following we study in more detail the dynamics of an isolated buoyant bubble or thermal, as it shall also be referred to.

## 3 Dynamics of isolated thermals

In the following we consider cases with different degrees of stratification and different extents of the computational domain.

### 3.1 Modest stratification

We begin by considering first vertical cross-sections of entropy and velocity; see Fig. 3. The archimedian buoyancy force is largest in the middle of the bubble, and that is also where the vertical velocity is largest. On both sides of the bubble the velocity turns over, as expected (compare with observations of thermals described by Scorer 1957).

A rather different impression is obtained when looking at the bubble in a comoving frame of reference; see Fig. 4. In this frame there is a stagnation point and hence there is a clear distinction of regions inside and outside the bubble. The flow pattern generally conforms with the notion of bubbles behaving like balloons with a more-or-less well defined surface and gas flowing around this surface. However, the bubble clearly grows in size and even its mass grows during its ascent.

It should be noted that this addition of new material with ‘different’ specific entropy does not contradict the bubble concept. This is because there is no entropy discontinuity between the bubble and the surroundings, so that the entropy of a captured particle can be changed gradually – by friction and by thermal diffusion – as it enters the bubble.

By the time $`t=30`$ the effects of the lateral boundaries have begun to affect the evolution of the bubble. Therefore we show in Fig. 5 the case of a wider box ($`L_x=8`$). Note that there is now a noticeable flow speed even beyond $`|x|=4`$ (the extent of the box in the previous case). In this calculation we have also included strong velocity perturbations, but the overall flow pattern is still dominated by the rising bubble.

In order to quantify the rise and the growth of the bubble in detail, we define the bubble $`B`$ as all points in space where $`ss_{\mathrm{crit}}`$, with $`s_{\mathrm{crit}}`$ just a little larger than the background value, which is here $`s_0=0.596`$, so we chose $`s_{\mathrm{crit}}=0.597`$. The volume of the bubble is then estimated as

$$V(t)=_B𝑑V,$$
(20)
and the volume radius of the bubble is

$$R(t)=\left(V/\frac{4}{3}\pi \right)^{1/3}.$$
(21)
The mass of the bubble is

$$M(t)=_B\rho 𝑑V.$$
(22)
In Fig. 6 we plot $`R(t)`$ and $`M(t)`$. Both functions increase monotonically, except for some minor departures at late times when the bubble has reached the top of the layer. Note also that at early times ($`t<4`$) $`R`$ increases somewhat faster than at later times. Qualitatively this type of behaviour is expected, because radiative diffusion causes structures to grow proportional to $`t^{1/2}`$, which causes an infinite slope of $`R(t)`$ at $`t=0`$.

In order to make further comparison with Hazlehurst’s theory of buoyant bubbles, we measure the position of the centre of mass of the bubble,

$$z_{\mathrm{bubble}}=_Bz\rho 𝑑V/_B\rho 𝑑V.$$
(23)
Since $`z`$ decreases upwards we define the height of the bubble as $`h=z_{\mathrm{bubble}}`$. The height $`h`$ and velocity $`v=\mathrm{d}h/\mathrm{d}t`$ of the bubble are shown in Fig. 7. Note that the height seems to approach a maximum near $`h=3.5`$, so the centre of mass of the bubble does note quite reach the top of the box. (Below we shall show that for larger bubbles, $`R_0=1`$, the maximum height is even less, suggesting that this is at least partly a geometrical effect; we shall also see that the top of the bubble does reach the top of the box in all cases.)

The momentum of the bubble, $`Mv`$, is plotted in Fig. 8 and compared with the time integrated buoyancy force,

$$_0^t(M_\mathrm{d}M)gdt,$$
(24)
where

$$M_\mathrm{d}(t)=_B\stackrel{~}{\rho }𝑑V$$
(25)
is the mass of displaced material and $`\stackrel{~}{\rho }`$ is the density of the undisturbed medium. Between $`t=5`$ and 15 the momentum of the bubble is somewhat larger than expected from the buoyancy force. This discrepancy depends somewhat on the definition of the boundary of the bubble. However, more dramatic is the sudden loss of momentum of the bubble after $`t20`$, whilst the buoyancy force, as estimated by Eq. (24), continues to operate beyond this time.

There are at least two possible reasons for this sudden braking effect. One reason could be that the blob gets too close to the top and loses momentum simply because of pressure build-up between the bubble and the top boundary. A second possibility could be that the bubble loses momentum due to some genuine resistance mechanism, such as wave braking or viscous friction. However, the braking effect seen in the simulations is too sudden and too strong to be explained by any genuine braking mechanism. Thus, we now turn to the first possibility, of which we can distinguish two variants. It is possible that the pressure build-up near the top is either an artifact of the top boundary being impenetrable, or it could be simply a feature of strong density stratification which causes the bubble to expand rapidly sideways. In order to drive strong sideways motions there must naturally be a horizontal pressure gradient which would also act in the vertical direction and slow down the ascent. This mechanism is known in compressible convection as buoyancy braking (Hurlburt et al. 1984).

In order to clarify the nature of the additional braking effect seen in the simulations, we first compare with a simulation using a somewhat taller box to see whether or not the braking sets in later, as would be the case if the impenetrable top boundary was the reason for the braking effect. (In the following we use calculations where $`R_0=1`$.)

### 3.2 Moving the top boundary further away

With $`\xi _0=0.3`$ and $`z_{\mathrm{top}}=6.0`$ we have $`z_{\mathrm{}}=6.75`$, see Eq. (8), so we can move the top boundary upwards by no more than about 10%. In the following, we discuss a model with $`z_{\mathrm{top}}=6.5`$, but otherwise the same stratification. This means that at the new boundary we have to change $`\xi _0`$ by $`\mathrm{\Delta }\xi =0.4\times \mathrm{\Delta }z=0.2`$, so we have to require $`\xi _0=0.1`$.

In Fig. 9 we compare the vertical entropy profiles along the axis of the bubble for $`z_{\mathrm{top}}=6.5`$ (with $`\xi _0=0.1`$) and $`z_{\mathrm{top}}=6.0`$ (with $`\xi _0=0.3`$). We also compare with the case of an open top boundary that we have modelled by putting an extra layer on top of the box where gravity goes smoothly to zero and, as in Brandenburg et al. (1996), radiative diffusion is replaced by a heating/cooling term of the form $`\tau ^1(z)\rho (TT_{\mathrm{top}})`$, where $`\tau ^1=0`$ everywhere except above $`z=z_{\mathrm{top}}`$ where it goes smoothly to $`\tau ^1=10`$. This procedure allows the flow to penetrate the layer $`z=z_{\mathrm{top}}`$ freely.

It turns out that the entropy profiles at the location of the bubble are not significantly affected by the properties of the top boundary. The entropy drop near the surface is a consequence of fixing the top temperature, $`T_{\mathrm{top}}`$. Any increase in the logarithmic pressure at the top, $`\delta \mathrm{ln}p_{\mathrm{top}}`$, causes a corresponding decrease in the entropy, $`\delta s_{\mathrm{top}}=0.4\delta \mathrm{ln}p_{\mathrm{top}}`$. Note, however, that the location of this entropy drop at the surface moves further away from the location of the bubble as we extend the box.

In Fig. 10 we show the pressure fluctuations (relative to the horizontal mean) together with velocity vectors. Near the top of the bubble there is a strong local maximum of the pressure fluctuation that drives the gas sideways. We have checked that the ram pressure integrated over the projected surface is roughly what is needed to explain the discrepancy between acceleration and buoyancy force. This is suggestive of buoyancy braking being the cause of the sudden drop of momentum seen in Fig. 8 (for $`R_0=0.5`$) and in Fig. 11 (for the present case of $`R_0=1`$).

In Fig. 11 we compare the evolution of $`h`$ in the two cases with different values of $`z_{\mathrm{top}}`$. Within the range of accuracy the two curves are consistent. Of course, the difference in the value of $`z_{\mathrm{top}}`$ is not very large, but the increase in the total number of scale heights covered in the simulation is significant: the value of $`\mathrm{ln}(p_{\mathrm{bot}}/p_{\mathrm{top}})`$ has increased from 6.2 to 9.0 pressure scale heights.

Note that $`h`$ reaches a maximum around $`3`$, i.e. significantly less than the value of $`|z_{\mathrm{top}}|`$. This is partly a geometrical effect, because smaller bubbles are able to travel somewhat further (in the previous subsection, where $`R_0=0.5`$ instead of 1.0, the bubble went to $`h=3.5`$). In the case of the open top boundary, the bubble rises till $`h=4`$, which is still small compared with $`|z_{\mathrm{top}}|=6`$. The relatively small values of $`\mathrm{max}(h)`$ are partly due to the fact that $`h`$ measures the position of the centre of mass, but in the strongly stratified case most of the mass is at the bottom of the bubble. Furthermore, in all cases considered the bubbles take mushroom form with a significant portion of the mass residing in the stem of the mushroom. A somewhat better representation of where most of the hot material resides is gained by looking at the value of the entropy weighted height, $`h_s=z_s`$, where

$$z_s=_Bz(ss_{\mathrm{crit}})\rho 𝑑V/_B(ss_{\mathrm{crit}})\rho 𝑑V,$$
(26)
which is also plotted in Fig. 11. However, even the entropy weighted height of the bubble is not very close to the top of the bubble. Nevertheless, the location $`h_{\mathrm{top}}`$ of the top of the bubble, i.e. where $`s=s_{\mathrm{crit}}`$, does reach value close to $`|z_{\mathrm{top}}|`$; see Fig. 11. There remains however some worry that the centre of mass of the bubble is generally unable to travel great distances. In order to clarify this possibility we now consider the case of weak stratification where we can easily increase the extent of the box.

### 3.3 Weak stratification

We now consider a case of weak stratification and two different values of $`z_{\mathrm{top}}`$ ($`6`$ and $`14`$); hence we choose $`\xi _0=30`$ and $`\xi _0=26.8`$ respectively, so that the background stratification is the same for $`6<z<2`$. The total number of pressure scale heights in these two cases is $`\mathrm{ln}(p_{\mathrm{bot}}/p_{\mathrm{top}})=0.25`$ and 0.5, respectively.

In Fig. 12 we compare the momentum balance for the two cases with different values of $`z_{\mathrm{top}}`$. Again, the agreement between the momentum and the time integrated buoyancy force is good up to $`t=10`$ (for $`z_{\mathrm{top}}=6`$) or up to $`t=30`$ (for $`z_{\mathrm{top}}=14`$). The presence of the boundary clearly influences the motion of the bubble; nevertheless the top of the bubble does manage to reach the boundary, as we see from Fig. 13 (for weak initial perturbations) and Fig. 14 (for strong initial perturbations). This is in contrast to the results of Sinjab et al. (1990) who conclude that the blobs in their calculations ‘fail to reach the top of the envelope’. In our calculation the hot part of the bubble seems to stop at some point just below the boundary and stays there. Furthermore, Sinjab et al. (1990) interpret their two-dimensional calculations as indicating a fragmentation. We found no indication of any such effect, irrespective of whether the stratification was weak or strong. For comparison we also carry out two-dimensional cartesian calculations verifying that in the two-dimensional case the eddies do eventually travel downwards when they come sufficiently close to the horizontal boundaries. Again, in the weakly stratified case the bubbles can travel to large heights, but the effects of the boundaries begin to become important much earlier than in the three-dimensional case; see Fig. 15.

The results of Sinjab et al. (1990) are of course for the stratified case. However, our point is that it is the restriction to two dimensions (in cartesian geometry) that causes boundary effects to become extremely pronounced. Since stratification itself can also act as an effective boundary, it was necessary to go to the weakly stratified case where it is possible to move the boundaries much further away. Although the box shown in Fig. 15 was big enough to prevent the edge of the bubble from moving down again, boundary effects did begin to affect the evolution at times as early as $`t=20`$.

The ability of our bubble to rise right to the top of the box might seem surprising when one recalls that atmospheric thermals only manage to reach a certain maximum height. However one must not forget that the buoyancy force, relative to gravity, is here orders of magnitude greater than in the case of a typical atmospheric thermal (see e.g. Scorer 1957). A better analogy might therefore be that of an air bubble rising through water and eventually reaching the surface.

## 4 Calculation of the entrainment parameter

The introduction of a parameter to describe the entrainment of fluid by a rising ‘cloud’ was proposed by Morton et al. (1956). They introduced an equation of the form

$$\dot{V}=4\pi R^2\alpha _\mathrm{v}v,$$
(27)
which they regarded as describing ‘conservation of volume’. We have here rewritten their Eq. (16) using a slight change of notation. The ‘constant’ $`\alpha _\mathrm{v}`$ will be referred to hereafter as the volume entrainment coefficient. The above equation can be simplified to

$$\dot{R}=\alpha _\mathrm{v}v,$$
(28)
where $`R`$ is the volume radius.

We note that Eq. (27) is taken over in the work of Turner (1963), with $`\alpha _\mathrm{v}`$ becoming Turner’s ‘alpha’.<sup>1</sup><sup>1</sup>1We note that although $`\alpha `$ has the same meaning in Morton et al. (1956) and Turner (1963) the characteristic velocities of the bubbles are defined differently, and this is ‘compensated’ by the inclusion of a $`k`$-factor in the entrainment equation of Morton et al. ($`k=\text{ratio}`$ of mean to axial velocity). This would not matter, except that when comparing with experiment, Morton et al. believe the observations relate to the axial velocity whereas Turner takes them as referring to the mean velocity. This means that when comparing experimental results it is really the $`\alpha k`$ of Morton et al. which should be compared with the $`\alpha `$ of Turner.

Now the concept of ‘conservation of volume’ lying behind Eq. (27) will be unfamiliar to many physicists. We therefore thought it to be worthwhile to concentrate instead on the accretion of mass rather than of volume and to write

$$\dot{M}=4\pi R^2\stackrel{~}{\rho }\alpha _\mathrm{m}v,$$
(29)
where $`\dot{M}`$ is the mass entrainment rate, $`\stackrel{~}{\rho }`$ some average of the surrounding material density near the bubble, and $`\alpha _\mathrm{m}`$ the mass entrainment coefficient. Although an average of $`\stackrel{~}{\rho }`$ over the bubble surface might be more appropriate, it is more convenient (and adequate for present purposes) to use a volume average. This gives

$$\dot{M}=\frac{3}{R}M_\mathrm{d}\alpha _\mathrm{m}v$$
(30)
where, as before, $`M_\mathrm{d}`$ is the mass displaced by the bubble.

We can now use Eqs (28) and (30) to calculate $`\alpha _\mathrm{v}`$ and $`\alpha _\mathrm{m}`$, all other quantities being already known. The results are shown in Figs 16 and 17.

It is interesting to note that for those parts of the curves not influenced by the special effects discussed previously, the quantity $`\alpha _\mathrm{v}`$ is in reasonable agreement with the experimental values of 0.25–0.34 (Morton et al. 1956, Turner 1963), with the agreement being better in the more strongly stratified case. For smaller initial bubble radii, $`R_0=0.5`$, $`\alpha _\mathrm{v}`$ attains noticeably larger values. Nevertheless, for both bubble radii investigated $`\alpha _\mathrm{m}`$ is roughly comparable ($`\alpha _\mathrm{m}0.25`$). Thus, the assumption of Turner (1963) that the value of $`R_0`$ does not influence ‘alpha’ can be verified only for $`\alpha _\mathrm{m}`$.

Better agreement between the ‘experimental’ and ‘calculated’ values cannot perhaps be expected since the Reynolds numbers for our calculations are lower than those appropriate to the experimental results, which leads to a more laminar entrainment process in our case.

Regarding $`\alpha _\mathrm{m}`$, we note that the entrainment of material leads to an effective drag force on the bubble given by

$$\text{Drag force}=\dot{M}v=4\pi R^2\stackrel{~}{\rho }\alpha _\mathrm{m}v^2,$$
(31)
which may be compared with the normal hydrodynamic drag

$$\text{Hydrodynamic drag}=\frac{1}{2}\pi R^2\stackrel{~}{\rho }C_\mathrm{D}v^2,$$
(32)
with $`C_\mathrm{D}`$ being the hydrodynamic drag coefficient. According to Moore (1967), normal hydrodynamic drag plays only a minor rôle in influencing the motion of a thermal. We can now check Moore’s point with the help of Eqs. (31) and (32); we confirm that hydrodynamic drag really does have only a minor influence on the motion, except possibly during the final stages in the weakly stratified case.

## 5 Details of the entrainment

Details of the entrainment process are viewed best in terms of tracer particles that are passively advected by the flow. In Fig. 18 we show the location of initially uniformly distributed particles at $`t=10`$. Particles that were originally inside the bubble (as defined by $`ss_{\mathrm{crit}}`$) remain inside the bubble for all times. However, an increasing portion of new particles from outside the bubble is constantly being entrained, which happens mostly through the top boundary of the bubble. These particles then move along the periphery of the bubble tailwards where they find their way into the centre of the ring vortex associated with the bubble.

In Fig. 18 we have also indicated the trajectory of three neighbouring particles that were originally above the bubble and were subsequently entrained and lifted upwards together with the bubble. The kinks in each of these trajectories correspond to the moment at which the particles were overtaken by the bubble and became entrained.

## 6 Conclusions

The main aim of this paper was to test the validity of the bubble or thermal concept; in this respect the following conclusions can be drawn.

We found that the bubble could easily be followed as a well-defined entity throughout the calculations. Nevertheless, its dynamical behaviour consisted of two distinct phases. In the first of these, the dynamical behaviour expected from the simple bubble dynamics of equating buoyancy forces to rate of momentum change was indeed confirmed. However, in the second phase an unexpected braking effect appeared. Further investigation led us to attribute this effect to a combination of boundary effects (artificial) and buoyancy braking (real).

We found, in contrast to the results of Sinjab et al. (1990), that the top part of the bubble goes on rising until it reaches the surface. Bubbles penetrating to the surface layers favour the view that dissipation is an important factor in understanding contact binaries (Hazlehurst 1985, 1996).

Since the code used by us was a three-dimensional one including effects of viscosity and thermal (radiative) conductivity, we were in a position to follow in detail and more realistically the mass changes of the bubble due to entrainment of material. This made it possible for us to calculate the value of the entrainment coefficient $`\alpha `$ rather than, as was previously necessary, to take some value from glider observations (atmospheric thermals) or experiment. The values found by us were in good agreement with the those derived empirically, although the non-constancy of ‘alpha’ means that we have throughout referred to the entrainment coefficient rather than the entrainment constant. Finally, the introduction of a mass entrainment coefficient to replace (or at least supplement) the previously used volume coefficient appears to us desirable on physical grounds.

###### Acknowledgements.

We thank Dr Robert Smith for useful comments regarding Fig. 18. We also thank an anonymous referee for suggestions which led to an improvement of the presentation. John Hazlehurst is grateful for the hospitality of Nordita during this investigation. We acknowledge the use of GRAND, a high-performance computing facility funded by PPARC.


format: text

text:
# Self-consistent solutions of Ginzburg–Landau equations and superconducting edge-suppressed states in magnetic field

## I Introduction

The behavior of finite size superconductors in magnetic field was studied on the base of nonlinear system of Ginzburg–Landau equations  in numerous theoretical papers (see, for instance, \[2–15\]). The results of \[2–15\] were used, in particular, to explain some anomalies, observed in a number of experiments with small size superconductors, placed in external magnetic field \[16–24\].

In our recent paper , by using numerical methods, the self-consistent solutions of nonlinear Ginzburg-Landau equations were studied, for a long superconducting cylinder of finite radius in external axial magnetic field, with no vortices inside the cylinder ($`m=0`$). It was found, that type-II superconducting cylinder, which is in the Meissner (vortex-free) state, upon increasing the field $`H`$, may pass by a jump into a special (also vortex-free) edge-suppressed state. The order parameter $`\psi `$ in this state is strongly suppressed (practically to zero) in the vicinity of the cylinder surface, however, the superconductivity survives near the specimen center, where $`\psi `$ remains finite. Such edge-suppressed solutions were not studied previously \[1–15\].

Apart from the the Meissner-state ($`m=0`$) more general one-dimensional solutions exist, when the vortex line, carrying arbitrary number of flux quanta ($`m=1,\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}3},\mathrm{}`$), is situated on the cylinder axis. Such solutions (which are called ”the giant-vortex states”) have been extensively studied in the literature (\[2-30\]) and used in the interpretation of experiments. The shape of the giant-vortex changes smoothly with the external field $`H`$. As is shown in the present paper, for a cylinder of finite radius, the giant-vortex state (with vorticity $`m>0`$) becomes unstable, when the field $`H`$ exceeds some critical value $`H_s^{(m)}`$. After that, the giant-vortex solution transforms its shape by a jump and acquires the ”edge-suppressed” form with the same vorticity $`m`$. Such edge-suppressed states are characterized by a strongly suppressed value of the order parameter $`\psi `$ near the cylinder boundary, and are analogous to the edge-suppressed (or, rim-suppressed) states, found in  for the case $`m=0`$. Previously, such edge-suppressed modifications of the giant-vortex states were not studied in detail, though some evidence of the edge-suppression was noticed also in .

We should stress, that the edge-suppressed state is noting else, but the usual giant $`m`$-vortex solution in a field $`H>H_s^{(m)}`$. If the field $`H`$ is increased further, the edge-suppressed solution continues to vary its form, the maximal value of the order parameter $`\psi _{max}^{(m)}`$ gradually diminishes, and for $`\psi _{max}^{(m)}1`$ the edge-suppressed solution may be expressed in terms of the Kummer functions (see, for instance, ).

As is shown below, the edge-suppressed states can be formed only in type-II superconductors, with $`\text{ϰ}>1/\sqrt{2}`$, where ϰ is the Ginzburg-Landau parameter. In type-I superconducting cylinder, with $`\text{ϰ}<1/\sqrt{2}`$, the edge-suppressed states do not exist.

It is known, that for $`m>0`$ the order parameter vanishes in the vicinity of the cylinder axis according to the law $`\psi (r)r^m`$ ; if the field increases, $`\psi `$ vanishes also near the surface (where the edge-suppressed state forms), so, the superconductivity persists only at some distance from the cylinder axis. If the magnetic quantum number $`m`$ is increased, the region, where the order parameter $`\psi (r)0`$, shifts more and more toward the cylinder surface. For a fixed cylinder radius ($`R=`$const, $`\text{ϰ}>1/\sqrt{2}`$) there exists some maximal value $`m=m_{max}`$, for which the edge-suppressed state degenerates into the usual state of the surface superconductivity \[26–30\].

The transition of a superconducting cylinder into the edge-suppressed state may be accompanied by jumps in the dependencies versus $`H`$ of such quantities, as: the specimen magnetic moment ($`4\pi M`$); the total magnetic flux, captured inside the system ($`\mathrm{\Phi }_1`$); the system Gibbs free energy ($`G`$); and by peculiarities in the behavior of other parameters, which characterize the superconducting state. These topics are discussed below in detail. No immediate comparison of the theory with the experiment is possible, because the model case of an infinitely long cylinder approximates the real samples geometry \[16-24\] very remotely. However, a number of qualitative predictions, which follow from this model, may be used for interpreting some of the peculiarities, observed in the experiments with mesoscopic samples. In particular, the controversial paramagnetic Meissner effect in superconductors is shortly discussed below, basing on the Ginzburg–Landau theory approach.

The paper is organized as follows. In Sec. 2, the basic equations and the boundary conditions of the problem are written. The numeric algorithm, used to find the self-consistent solutions of nonlinear system of equations, is shortly described. Sec. 3 contains the results of numerical calculations. We are forced to present our results by a large number of graphics, which illustrate different details of the system behavior in this many parameter problem (the solutions depend on $`m`$, ϰ, $`R`$, $`H`$, the specimen temperature $`T`$, the critical temperature $`T_c`$ and the coherence length $`\xi _0`$). The necessary comments accompany the presentation of the material. Sec. 4 contains a short resume and discussion of the results.

## II The setting of the problem

Consider long superconducting cylinder of radius $`R`$ in the external magnetic field $`H`$, which is parallel to the cylinder element. The basic system of the Ginzburg-Landau equations  is of the form

$$\mathrm{rot}\mathrm{rot}𝐀=\frac{4\pi }{c}𝐣_s,\frac{4\pi }{c}𝐣_s=\frac{\psi ^2}{\lambda ^2}\left(\frac{\varphi _0}{2\pi }\mathrm{\Theta }𝐀\right),$$

$$^2\psi \left(\mathrm{\Theta }\frac{2\pi }{\varphi _0}𝐀\right)^2\psi +\frac{1}{\xi ^2}(\psi \psi ^3)=0,$$

where $`𝐀`$ is the vector-potential of the magnetic field ($`𝐁=\mathrm{rot}𝐀`$), $`𝐣_s`$ is the current density inside the superconductor, $`\lambda `$ is the field penetration depth, $`\xi `$ is the coherence length, $`\lambda =\text{ϰ}\xi `$. The order parameter, in a general case, is written as $`\mathrm{\Psi }=\psi e^{i\mathrm{\Theta }}`$, where $`\psi `$ is the modulus and $`\mathrm{\Theta }`$ is the phase of the order parameter. From the single-valuedness of $`\mathrm{\Psi }`$ the condition follows

$$_C\mathrm{\Theta }d𝐥=2\pi m,$$

where the contour $`C`$ embraces the vortex axis, $`m`$ is integer (the topological invariant, or fluxoid, or vorticity), which shows, how many vortices may be present inside the contour $`R`$.

In the cylidrical system of co-ordinates $`r,\phi ,z`$, with $`z`$ axis directed along the cylinder element \[when the vector-potential has only one component, $`𝐀=𝐞_\phi A(r)`$\], these equations may be written in the dimensionless form

$$\frac{d^2U}{d\rho ^2}\frac{1}{\rho }\frac{dU}{d\rho }\psi ^2U=0,$$
$`(1)`$
$$\frac{d^2\psi }{d\rho ^2}+\frac{1}{\rho }\frac{d\psi }{d\rho }+\text{ϰ}^2(\psi \psi ^3)\frac{U^2}{\rho ^2}\psi =0.$$
$`(2)`$
Here, instead of the dimensioned potential $`A`$, field $`B`$ and current $`j_s`$, the dimensionless quantities $`U(\rho ),b(\rho )`$ and $`j(\rho )`$ are introduced:

$$A=\frac{\varphi _0}{2\pi \lambda }\frac{U+m}{\rho },B=\frac{\varphi _0}{2\pi \lambda ^2}b,b=\frac{1}{\rho }\frac{dU}{d\rho },$$
$`(3)`$
$$j(\rho )=j_s/\frac{c\varphi _0}{8\pi ^2\lambda ^3}=\psi ^2\frac{U}{\rho },\rho =\frac{r}{\lambda }.$$

(The field $`B`$ in (3) is normalized by $`H_\lambda =\varphi _0/(2\pi \lambda ^2)`$, with $`b=B/H_\lambda `$; instead of $`H_\lambda `$ one can normalize by $`H_\xi =\varphi _0/(2\pi \xi ^2)H_{c2}`$ \[28–30\], or by thermodinamical critical field $`H_c=\text{ϰ}H_\lambda /\sqrt{2}`$. The coefficents in (1)–(3) would change accordingly.)

The magnetic flux, confined inside the contour of the radius $`r`$, is

$$\mathrm{\Phi }=𝐁𝑑𝐬=_r𝐀𝑑𝐥=\varphi _0(U+m),U=U(\rho ),\rho =\frac{r}{\lambda }.$$

Thus, the potential $`U(\rho )`$, in the normalization adopted above, is related to the flux $`\mathrm{\Phi }(\rho )`$ by a simple formula $`\varphi \mathrm{\Phi }/\varphi _0=U(\rho )+m`$.

Because the magnetic flux through the vanishing contour is zero, and the field $`B|_{r=R}=H`$, the following boundary conditions correspond to Eq. (1):

$$U|_{\rho =0}=m,\frac{dU}{d\rho }|_{\rho =\rho _1}=h_\lambda ,$$
$`(4)`$
where $`\rho _1=R/\lambda `$, $`h_\lambda =H/H_\lambda `$, $`H_\lambda =\varphi _0/(2\pi \lambda ^2)`$.

As to the Eq. (2), we shall take the usual boundary condition on the external surface : $`d\psi /d\rho |_{\rho =\rho _1}=0`$. The order parameter at the center is either maximal (if $`m=0`$), or zero (if $`m>0`$, see, for instance, ), thus, the following boundary conditions correspond to Eq. (2):

$$\frac{d\psi }{d\rho }|_{\rho =0}=0,\frac{d\psi }{d\rho }|_{\rho =\rho _1}=0(m=0),$$
$`(5)`$
$$\psi |_{\rho =0}=0,\frac{d\psi }{d\rho }|_{\rho =\rho _1}=0(m>0).$$
$`(6)`$
The magnetic moment (or, magnetization) of the cylinder, related to the unity volume, is

$$\frac{M}{V}=\frac{1}{V}\frac{BH}{4\pi }𝑑v=\frac{B_{av}H}{4\pi },$$
$`(7)`$
$$B_{av}=\frac{1}{V}B(𝐫)𝑑v=\frac{1}{S}\mathrm{\Phi }_1,$$

where $`B_{av}`$ is the mean field value inside the superconductor, $`\mathrm{\Phi }_1=\mathrm{\Phi }(R)`$, $`S=\pi R^2`$. In normalization (3), denoting $`\overline{b}=B_{av}/H_\lambda `$, $`h_\lambda =H/H_\lambda `$, $`M_\lambda =M/H_\lambda `$, one finds from (7):

$`4\pi M_\lambda =\overline{b}h_\lambda ,\overline{b}={\displaystyle \frac{2}{\rho _1^2}}(U_1+m),\varphi _1={\displaystyle \frac{\mathrm{\Phi }_1}{\varphi _0}}=U_1+m,`$
$`U_1=U(\rho _1),\rho _1={\displaystyle \frac{R}{\lambda }}.(8)`$
For the difference of Gibbs free energies of the system in superconducting and normal states, $`\mathrm{\Delta }G=G_sG_n`$, it is convenient to use the exact expression:

$$\mathrm{\Delta }G=_{s0}\frac{1}{2}MH+\frac{\varphi _0Z}{8\pi }m(B(0)H),$$
$`(9)`$
$$_{s0}=\frac{H_c^2}{8\pi }\left[\psi ^42\psi ^2+\xi ^2\left(\frac{d\psi }{dr}\right)^2\right]𝑑v,H_c=\frac{\varphi _0}{2\sqrt{2}\pi \lambda \xi },$$

where $`_{s0}`$ corresponds to the superconductor condensation energy, $`M`$ is the cylinder magnetic moment, $`Z`$ is the cylinder length, $`m`$ is the magnetic quantum number, $`B(0)`$ is the magnetic field at the cylinder axis, $`H`$ is the external magnetic field. The expression (9) follows from the general formula, obtained in  for the free energy of a hollow superconducting cylinder, if the radius of a hollow is put to zero (see also ). Using normalization (3), and also (8), one finds from (9) the normalized expression

$$\mathrm{\Delta }g=\mathrm{\Delta }G/\left(\frac{H_{c\mathrm{m}}^2}{8\pi }V\right)=g_0\frac{8\pi M_\lambda }{\text{ϰ}^2}h_\lambda +\frac{4m}{\text{ϰ}^2}\frac{b(0)h_\lambda }{\rho _1^2},$$
$`(10)`$
$$g_0=\frac{2}{\rho _1^2}_0^{\rho _1}\rho 𝑑\rho \left[\psi ^42\psi ^2+\frac{1}{\text{ϰ}^2}\left(\frac{d\psi }{d\rho }\right)^2\right].$$

The expressions (3),(8),(10) will be used in Sec. 3 for calculating the corresponding quantities.

We remind, that the field penetration depth $`\lambda `$ and the coherence length $`\xi =\lambda /\text{ϰ}`$ depend on temperature. Thus, the expressions above depend implicitly on temperature, and, formally, are valid for arbitrary values of $`T`$. (Though, the Ginzburg-Landau equations themselves are applicable only in the limit $`TT_c`$ \[34,28–30\], when the expression $`\xi (T)=\xi _0/\sqrt{1T/T_c}`$ may be used for the coherence length).

It is appropriate to make here a comment on the iteration procedure, we used to obtain the self-consistent solution of the system of equations (1)–(6). First, some trial function $`\psi (\rho )`$ was chosen, and the solution for $`U(\rho )`$ was found from Eqs. (1), (4). Then, this function $`U(\rho )`$ was introduced into Eq. (2), which was solved with account of the boundary conditions (5), (6) and new function $`\psi (\rho )`$ was found. Further, Eq.(1) was solved again, and all the procedure was repeted, until the functions $`\psi (\rho )`$ and $`U(\rho )`$ ceased to change, representing the self-consistent solution of the system (see also ). However, the numeric algorithm, used in our calculations, does not allow to find the unstable ”superheated” and ”supercooled” states in type-I superconductors (with $`\text{ϰ}<1/\sqrt{2}`$), where other methods are required (see ).

The question of stability requires an additional illucidation. Alongside with the centrally symmetric giant-vortex states, the so-called multi-vortex solutions of the same vorticity $`m`$ exist. The multi-vortex states correspond to the case, when a single giant vortex (with vorticity $`m`$) decays into a number of vortices (with vorticities $`m_i,m_i=m`$). These vortices can be arbitrary positioned on the cylinder cross-section, so the multi-vortex state may be no longer centrally symmetric. The question, which of these multi-vortex states possesses the minimal free energy, is important, because the physical system is usually assumed to occupy the state with the smallest energy (the equilibrium ground state), and all other states are considered as energetically unstable (or metastable). The free energy calculations to find the equilibrium multi-vortex configuration is a difficult problem, which requires the specific methods of investigation. (See, for instance, the theory , where the triangular lattice of single vortices is found as an equilibrium state for the case of infinite superconductor, placed in the field $`HH_{c1}`$. The numeric calculations for mesoscopic suprconducting discs may be found in \[6–15\], where a simplified model of the vortex distribution is considered. To find the equilibrium vortex configuration in a general case of a finite-size superconductor, placed in arbitrary magnetic field $`H`$, is a problem as yet unsolved.)

In the present work a limited task was set: to describe the formal properties of the giant-vortex states, which were not reported previously (in particular, the edge-suppressed transitions and the jumps in the magnetization curves with a fixed vorticity $`m`$). The general problem of energetic stability and of the equilibrium vortex configuration is left outside the scope of the present investigation. Some of the edge-suppressed states, which are described below, might in reality be energetically metastable and decay into an equilibrium ground state configuration. However, as is well known, the physical system not necessarily must be in equilibrium, but (with some probability) may also occupy the excited metastable states. The giant-vortex and the edge-suppressed solutions provide the system this opportunity, so the behavior of a superconductor in such states might be of some interest.

The results of our calculations are presented below.

## III The numerical results

The results of the investigation are illustrated below by a number of graphics, which detalize various aspects of this multy-parameter task. (Some of these results can be obtained by the methods, used in the papers, cited above. However, for the uniformity of presentation, we prefer to reproduce them by our method.)

Evidently, the superconductivity would be destroyed by a large external field $`H`$ with subsequent transition to the normal state. In Fig. 1($`a`$) the phase diagram is depicted (for $`\text{ϰ}=2`$ and different $`m`$), which separates the normal ($`n`$) and superconducting ($`s`$) states on the plane of variables $`R_\lambda =R/\lambda `$ and $`h_\xi =H/H_\xi `$. The phase diagram is a many-sheet surface, each sheet corresponds to the fixed value of the vorticity $`m`$. The sheet $`m=0`$ is the lowest laying (according to the increasing $`m`$), the sheet $`m=1`$ lays above it, still more above lays the sheet $`m=2`$, and so on. (It is sufficient to consider only the values $`m0`$, because the results for $`m<0`$ are easily reproducible with account of the evident symmetry $`(m,H)(m,H)`$.)

On each of $`m`$-sheets of the phase surface there is a phase curve ($`h_c`$, which may be marked by the index $`m`$), which separates the region of values $`(R_\lambda ,h_\xi `$), where the superconducting solutions exist (i.e. $`\psi >0`$), from the region of normal state solutions ($`\psi 0`$). For a bulk superconductor ($`R_\lambda 1`$) the superconducting state boundary (in the case of type-II superconductor) lays at $`h_\xi =H/H_\xi =1`$, i.e. at $`H=H_{c2}=\varphi _0/(2\pi \xi ^2)`$ \[28–30\] ($`H_{c2}H_\xi `$). As can be seen from Fig. 1($`a`$), for small $`R_\lambda `$ the superconducting states (with different values of $`m`$) are also possible in the fields, exceeding the critical value $`H_{c2}`$ (i.e. for $`h_\xi >1`$). For fixed $`R_\lambda `$ and $`m`$ there exists a maximal field $`h_\xi =h_c^{(m)}`$, above which the superconducting solution is impossible. As is well known, the maximal field, at which the superconductivity is still possible (for a bulk specimen with $`R_\lambda 1`$) is $`h_\xi =1.69`$ (i.e. $`H=H_{c3}`$, where $`H_{c3}=1.69H_{c2}`$ is the field, at which the surface superconduvtivity nucleates in macroscopic specimens \[26–30\]). As is evident from Fig. 1($`a`$), for small $`R_\lambda `$, the superconductivity is also possible in the fields, exceeding $`h_\xi =1.69`$. \[The analogous phase diagram in the case of type-I superconductor (with $`\text{ϰ}=0.5`$, $`m=0,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2}`$) is depicted in Fig. 1($`b`$).\]

Inside each of the superconducting regions in Fig. 1($`a`$) there exists a singular line $`h_s`$, which is depicted, for clarity, on a separate phase $`m`$-sheet (see Fig. 2). In Fig. 2($`a`$) the phase sheet $`m=0`$ (for $`\text{ϰ}=2`$) is depicted. The solid line $`h_c`$ corresponds to the superconducting state boundary. The letter $`s`$ marks the region, where the superconducting state exists ($`\psi 0`$), the letter $`n`$ denotes the normal state ($`\psi =0`$).

To every point($`R_\lambda ,h_\xi `$), which lays inside the superconducting region in Fig. 2($`a`$), corresponds some solution of Eqs. (1)–(6). When this point is shifted, the solution changes accordingly. The central part of the drawing (which lays in the withinity of $`h_\xi 0`$) corresponds to the superconducting Meissner state. The order parameter here is $`\psi 1`$, the external field $`H`$ is almost completely screened out and penetrates only to the distance $`\lambda `$ from the superconductor surface. The shape of the Meissner solution depends on $`H`$ only slightly.

There exists a critical line $`h_s`$ in Fig. 2($`a`$), near which the Meissner solution ($`\psi 1`$) starts to change strongly. The character of this transformation depends on the cylinder radius $`R_\lambda `$. If the representation point $`(R_\lambda ,h_\xi )`$ crosses the the critical line $`h_s`$ above the point $`R_\mathrm{\Delta }`$ ($`R_\lambda >R_\mathrm{\Delta }`$), the static Meissner solution becomes absolutely unstable relative small variations in its form, and transforms by a jump into new stable (static) superconducting solution with the same $`m`$. In this new state (which may be called as ”edge-suppressed”, or ”rim-suppressed” ) the order parameter $`\psi (r)`$ is strongly suppressed in some layer, which is situated near the superconductor surface, while the magnetic field penetrates this layer practically without screening (see Fig. 2 in  for $`m=0`$ and also Figs. 3–5 below).

The jump transformation of the order parameter shape is acompanied by jumps in the magnetization and in other physical characteristics (see Fig. 3 below). The amplitude of these jumps diminishes with $`R_\lambda `$ diminishing, and for $`R_\lambda =R_\mathrm{\Delta }`$ the jump amplitude (see the broken vertical lines in Fig. 3) vanishes. Thus, for $`R_\lambda >R_\mathrm{\Delta }`$ the transition of the Meissner state ($`m=0`$) to the edge-suppressed form ($`m=0`$) is the first order phase transition (by a jump).

For $`R_\lambda <R_\mathrm{\Delta }`$ the Meissner solution also changes its form significantly, when the line $`h_s`$ is crossed, but this transformation happens smoothly, without jump (a second order phase transition to the edge-suppressed form). In the case $`R_\lambda =R_\mathrm{\Delta }`$ the jump amplitude is zero, so the first and second order phase transitions to the edge-suppressed state become indistinguishable. (The meaning of the points $`R_w,R_0,R_\mathrm{\Delta }`$ in Figs. 2($`a,d`$) will be explained later.)

As was noted, to every point of $`s`$-region in Fig.2($`a`$) corresponds some solution of Eqs. (1)–(6) for functions $`\psi `$ and $`U`$, which describe the superconducting state. When the representation point ($`R_\lambda ,h_\xi `$) moves, this state changes. Let us trace in more details, what happens, if the parameter $`R_\lambda =R/\lambda `$ remains constant (i.e. the temperature remains constant), but the parameter $`h_\xi =H/H_\xi `$ changes (i.e. only the field $`H`$ changes).

In Fig. 3 are depicted as functions of $`H`$: ($`a`$) – the system magnetic moment ($`4\pi M_\lambda `$) (8); ($`b`$) – the average magnetic field $`\overline{b}`$ \[or, the total magnetic flux inside the cylinder, $`\varphi _1=\frac{1}{2}\rho _1^2\overline{b}`$ (8)\]; ($`c`$) – the difference of the system free energies, $`\mathrm{\Delta }g`$, (10); ($`d`$) – the magnetic field magnitude at the cylinder axis, $`b_0`$; ($`e`$) – the maximal value of the superconducting order parameter, $`\psi _{max}`$; ($`f`$) – the order parameter value at the cylinder surface, $`\psi _1`$. Fig. 3 corresponds to the case $`m=0`$, $`\text{ϰ}=2`$, $`R_\lambda =5`$. \[The superconducting region in Fig. 2($`a`$) is crossed along the line $`R_\lambda =5`$.\] The values of $`h_\lambda =\text{ϰ}^2h_\xi `$ are plotted on the horizontal axis. By open circles the points are marked, which deserve a special commentary.

As is seen from Fig. 3($`a`$), the magnetic moment ($`4\pi M_\lambda `$) increases linearly for small $`h_\lambda >0`$ ($`h_\lambda =0`$ at the point 4). The average magnetic field inside the cylinder, $`\overline{b}`$, is essentially reduced, in comparison to the normal state (see the dotted line $`n`$ in Fig. 3($`b`$)). (Analogously behaves the total magnetic flux, confined inside the cylinder, $`\varphi _1`$). The initial linear part reflects the Meissner effect, it corresponds to the external field expulsion from the superconductor interior. If the field $`H`$ is increased further, the magnetic moment increases up to the point 5, where the magnetization diminishes by jump, and takes the value 6.

The analogous jumps are present in other quantities, depicted in Figs. 3($`bf`$). The value of the order parameter at the surface, $`\psi _1`$ (Fig. 3($`f`$)), changes especially strongly: from $`\psi _1=0.706`$ at the point 5 (where $`h_\lambda =1.6837`$), to $`\psi _1=0.087`$ at the point 6 (where $`h_\lambda =1.6838`$). The jump in the free energy $`\mathrm{\Delta }g`$ (Fig. 3($`c`$)) is not seen in this scale.

When the external field $`h_\xi `$ is increased beyond the point 6, the maximal value of the order parameter, $`\psi _{max}`$, starts diminish quickly. This is acompanied by a fast drop in the ”tail” of the magnetization. (This ”tail” lays between the points 6 and 7 in Fig. 3($`a`$) and represents the edge-suppressed solution, $`e_0`$.) When the field $`h_\xi `$ reaches the point 7, the order parameter vanishes completely, by the second order phase transition to the normal state ($`\psi _{max}0`$).

Fig. 4 ($`m=0`$, $`\text{ϰ}=2`$) demonstrates, what is going on, when the cylinder radius $`R`$ is diminishing. \[The superconducting region in Fig. 2($`a`$) is crossed along the line $`R_\lambda =3`$.\] Now, there are no jumps on the magnetization curve (Fig. 4($`a`$)), which were represented by dashed lines in Fig. 3($`a`$). But, the points of inflection 2 and 4 remain, where the magnetization curve has the maximal derivative. The corresponding points in Figs. 4($`b,c`$) are marked by the same numerals 2 and 4 as in Fig. 4($`a`$). Evidently, the solution transforms smoothly, when the field $`h_\lambda `$ passes the critical value $`h_s`$ (the point 4). This point corresponds approximately to the value $`R_\mathrm{\Delta }`$ in Fig. 2($`a`$), where the derivative of the magnetzation curve is infinite.

In Fig. 5 the case $`m=0`$, $`\text{ϰ}=2`$, $`R_\lambda =1`$ is presented. \[The superconducting region in Fig. 2($`a`$) is crossed along the line $`R_\lambda =1`$.\] Here, the inflection point 4 ($`h_\lambda =5.91`$) on the curve ($`4\pi M_\lambda `$) practically coincides with the termination point 5 ($`h_\lambda =5.96`$). This situation corresponds to the point $`R_w`$ in Fig. 2($`a`$) (where $`h_s`$ and $`h_c`$ merge), and where the tail of the magnetization disappears (i.e. the width of the tail vanish). For $`R_\lambda <R_w`$ the second order phase transition to the normal state occures, and the magnetization curve has no tail.

Evidently, the picture, presented by Fig. 2($`a`$) ($`m=0`$), is symmetric relative the axis $`h_\xi =0`$. The analogous picture, which corresponds to the phase sheet $`m=1`$ \[Fig. 2($`d`$)\], is asymmetric relative the axis $`h_\xi =0`$. The edge-suppressed state exists only in the region $`h_\xi >0`$, and they do not exist in the region $`h_\xi <0`$. This may be interpreted in the physical terms.

The phase sheet $`m=1`$ describes situation, when at the cylinder axis there is a single vortex (its own field is assumed to have a positive sign), with the order parameter vanishing on the axis: $`\psi _0=0`$. If the external field of the same direction ($`h_\xi >0`$) is applyed, the magnetic field is additionally pumped into the specimen. When the field reaches the critical value $`h_s`$, the superconductivity is partly destroyed, and a new superconducting state forms, with strongly suppressed order parameter $`\psi _1`$ at the cylinder surface (see Fig. 6($`f`$) below). This state preserves up to the field $`h_\xi =h_c^{(+)}`$ (the right branch of the critical line, $`h_c>0`$), until the order parameter vanishes finally everywhere. Thus, in the vicinity of the line $`h_c^{(+)}`$ in Fig. 2($`d`$), the mechanism of the order parameter suppression by the external field acts.

If the external field has the opposite sign \[$`h_\xi <0`$, Fig. 2($`d`$)\], the field of the vortex ($`m=1`$) is gradual pumped out from the vortex interior. When the external field reaches the critical value $`h_\xi =h_c^{()}`$ (the left branch of the critical line, $`h_c<0`$), the vortex magnetic field $`b_0`$ is completely sucked out from the superconductor (see Fig. 6($`d`$)). When this happens, the vortex field at the axis vanishes (this means the jump transition from the state $`m=1`$ to the state $`m=0`$); the order parameter also vanishes by a jump (see Figs. 6($`e,f`$)). Thus, in the vicinity of the line $`h_c^{()}`$, the mechanism of pumping out the field from the vortex interior acts. This explains the asymmetric behavior of the curves $`h_c^{()}`$ and $`h_c^{(+)}`$ in Fig. 2($`d`$).

Note also, that the curves $`h_c^{()}`$ and $`h_c^{(+)}`$ come up to the axis $`h_\xi =0`$, having different derivatives (the fracture point at $`h_\xi =0`$ is marked by an open circle in Fig. 2($`d`$); the corresponding value at $`R_\lambda `$-axis is $`R_0`$).

Mention one more peculiarity of the curves, depicted in Fig. 2($`d`$). In the case of small radii (for $`R_\lambda <R_0`$ and $`m=1`$) there exists a field interval, where the superconducting state, which was impossible at smaller $`h_\xi `$, becomes possible again for larger $`h_\xi `$. Here one has an example of the so-called ”reentry” superconductivity. In a cylinder of very small radius ($`R_\lambda <R_0`$) the magnetic field of the vortex can not be confined inside the superconductor and dissipates outside through the cylinder surface. However, the imposition of the external field prevents the dissipation of the vortex field and stabilizes the superconducting state. If the external field is increased further, the superconductivity will be finally destroyed. Thus, on the lowest part of the curve $`h_c`$ in Fig. 2($`d`$) there exists a minimum at $`R_\lambda =R_{\lambda min}`$ (marked by an open circle).

If the parameter ϰ is increased, the curves $`h_s`$ and $`h_c`$ behave analogously to what is shown in Figs. 2($`a,d`$). But, if ϰ diminishes (see Figs. 2($`b,e`$)), the point, where the curves $`h_s`$ and $`h_c`$ merge, raises up to larger $`R_\lambda `$, and the width of the region, where the edge-suppressed state exists, diminishes. At $`\text{ϰ}=1/\sqrt{2}`$ the width of this region vanishes, and the curves $`h_s`$ and $`h_c`$ merge into one indivisible curve $`h_c`$. For $`\text{ϰ}<1/\sqrt{2}`$ (i.e. for type-I superconductors) the edge-suppressed states do not exist. The behavior of the critical curves $`h_c`$ in the case of type-I superconductor is depicted in Fig. 2($`c,f`$) for different $`\text{ϰ}<1/\sqrt{2}`$ ($`m=0`$ and $`1`$). One can see, that the superconducting state with a vortex ($`m=1`$) at the cylinder axis, in principle, can exist in both type-I and type-II superconductors.

The solution behavior in a case $`m=1`$, $`\text{ϰ}=2`$ \[when the superconducting region $`s`$ in Fig. 2($`d`$) is crossed along the line $`R_\lambda =5`$\] is illustrated in Fig. 6. There are also jumps, as in Fig. 3, but the curves are asymmetric relative the axis $`h_\xi =0`$. In particular, on the magnetic moment curve (Fig. 6($`a`$)), in a case $`h_\lambda <0`$ there is no ”tail”, which exists in Fig. 3($`a`$) for $`h_\lambda <0`$ (between points 1 and 2). This asymmetry (as was already mentioned) may be explained by different mechanisms, which act in destruction of the superconducting state $`m=1`$. If $`H>0`$, the additional pumping of the external field into the superconductor occures, with the subsequent suppresion of the order parameter by the field, followed by the transition to the normal state (or, to the state $`mm+1`$). If $`H<0`$, the field is pumped out from the vortex and the superconductor passes into the vortex-free state $`m=0`$ (or, even into the normal state).

The analogous picture (for $`m=1`$, $`\text{ϰ}=2`$, $`R_\lambda =3)`$ is presented in Fig. 7. In difference to Fig. 4 ($`m=0`$), the curves in Fig. 7 ($`m=1`$) have no smooth ”tails” for $`h_\xi <0`$ (to the left of points 2), but they terminate by a jump at the points 2, where the field at the vortex axis, $`b_0`$, vanishes, and the order parameter suffers a jump.

Fig. 8 corresponds to the case $`m=1`$, $`\text{ϰ}=2`$, $`R_\lambda =1`$. Notice, that in Fig. 8($`a`$), at positive values $`h_\lambda >0`$ (in the interval between the points 3 and 4) the magnetic moment is positive ($`4\pi M_\lambda >0`$), which corresponds to the paramagnetic susceptibility. In the interval between the points 5 and 6 the magnetic moment is negative ($`4\pi M_\lambda <0`$), which corresponds to the diamagnetic susceptibility. As can be seen from Fig. 8($`b`$), in the interval between points 3 and 4, the external field enhances the superconductivity (because $`\psi _{max}`$ increases). In the interval between points 4 and 5, the external field suppresses the superconductivity (because $`\psi _{max}`$ decreases).

The analogous dependencies are presented in Fig. 9 for type-I superconductor ($`R_\lambda =5`$, $`\text{ϰ}=0.5`$, $`m=0,\mathrm{\hspace{0.17em}1}`$). The singular points are marked by an open circles. There are no ”tails” on the magnetization curves (Fig. 9($`a`$)) in both regions $`h_\xi <0`$ and $`h_\xi >0`$. This is characteristic for type-I superconductors with $`\text{ϰ}<1/\sqrt{2}`$, when the destruction of the superconductivity by the magnetic field proceeds in a jump, as the first order phase transition. \[This is valid, if $`R_\lambda >R_\mathrm{\Delta }`$ in Fig. 1($`b`$). If $`R_\lambda <R_\mathrm{\Delta }`$, the phase transition to the normal state is of second order, even in type-I superconductors; see Fig. 14 below.\]

There are other features, common to all type-I superconductors. Thus, if the family of curves, depicted in Fig. 2($`e`$) (for $`m=0`$ and different $`\text{ϰ}<1/\sqrt{2}`$), is re-drawn on the ($`R_\xi ,h_\xi `$)-plain, one obtains the unified phase curve, Fig. 10($`a`$) ($`m=0`$). Notice also, that in the case $`m=0`$ the curves with $`\text{ϰ}<1/\sqrt{2}`$ and $`\text{ϰ}1/\sqrt{2}`$ almost coincide, so the phase curves in Fig. 10($`a`$) look like one universal function. (For $`h_\xi >1`$ this function is $`R_\xi 2.8/h_\xi `$).

Analogously, if the family of curves in Fig. 2($`f`$) is re-drawn on the ($`R_\xi ,h_\xi `$)-plain, the Fig. 10($`b`$) emerges. One can see from this figure, that for all $`\text{ϰ}1/\sqrt{2}`$ there exists a minimal radius $`R_{\xi min}1.3`$. For $`R_\xi <R_{\xi min}`$ the superconducting state with a vortex ($`m=1`$) on the cylinder axis is impossible.

In Fig. 10($`b`$) ($`m=1`$), among others, the phase diagrams are depicted for $`\text{ϰ}1`$. The corresponding curves (for $`h_\xi >0`$ and $`R_\xi 1`$) have the asymptote $`h_\xi =1`$ (or $`H=H_{c2}`$), common to all type-II superconductors with $`\text{ϰ}1`$ (see Figs. 1 and 2). However, if $`\text{ϰ}<1`$, these asymptotes depend on ϰ and began deviate from the value $`h_\xi =1`$. \[Thus, $`h_\xi 1`$ for $`\text{ϰ}=1`$; $`h_\xi 1.058`$ for $`\text{ϰ}=0.9`$; $`h_\xi 1.218`$ for $`\text{ϰ}=0.8`$; $`h_\xi 1.416`$ for $`\text{ϰ}=1/\sqrt{2}`$.\] This means, that in the case $`m=1`$ the value $`\text{ϰ}=1`$ is singular for equations (1)–(6). At this value the solution passes from one branch ($`\text{ϰ}1`$) to another ($`\text{ϰ}<1`$). At the same time, the value $`\text{ϰ}=1/\sqrt{2}`$ is also singular, because at $`\text{ϰ}=1/\sqrt{2}`$ the second branch of solutions appears (or vanishes), which describe the edge-suppressed states \[compare the behavior of the curves $`h_s`$ in Figs. 2($`ae`$)\].

Now, we present the examples of concrete solutions of Eqs. (1)–(6), as functions of the space coordinate $`\rho =r/\lambda `$ ($`0rR`$) for $`m=1`$. (The solutions profiles for the vortex-free state ($`m=0`$) are shown in , Fig. 2.)

In Fig. 11 the solutions are shown in the case $`R_\lambda =5`$, $`m=1`$, for several values of $`h_\lambda `$. Note the curve 3 in Fig. 11($`a`$), which corresponds to $`h_\lambda =0`$ and describes a screened Meissner-type vortex state $`v_1`$ ($`m=1`$). (In Fig. 6($`a`$) to this state corresponds the point 4.) The curve 4 in Fig. 11($`a`$) corresponds to the point 5 in Fig. 6($`a`$), where $`h_\lambda =1.8254`$. Here the solution branch $`v_1`$ terminates, the temporal instability develops, and at $`h_\lambda =1.8255`$ the solution passes to new stable (static) branch $`e_1`$ (the edge-suppressed state $`m=1`$, the point 6 in Fig. 6($`a`$)). The curve 5 in Fig. 11($`a`$) corresponds to the value $`h_\lambda =1.8255`$ and describes the edge-suppressed state $`e_1`$ with a vortex ($`m=1`$) at the cylinder axis (here $`\psi (0)=0`$ and the order parameter is additionally strongly suppressed near the cylinder surface).

It is expedient to trace in more details the system behavior in the states with different $`m`$, when the cylinder radius $`R`$ is fixed, and only the external field $`H`$ varies. If the superconducting region in Fig. 1($`a`$) is crossed along the line $`R_\lambda =2`$, the states with $`m=0,1,2,\mathrm{},11`$ are possible. These states are represented in Fig. 12 by the following dependencies: ($`a`$) – the normalized free-energy difference ($`\mathrm{\Delta }g`$); ($`b`$) – the magnetization ($`4\pi M_\lambda `$); (c) – the maximal value of the order parameter ($`\psi _{max}`$) in the state $`m`$.

In Fig. 12($`a`$) the free-energy curves with the adjacent $`m`$-values intersect at the points, marked by the open circles. At these points the equilibrium transitions from the state $`m`$ to the lower energy state $`m+1`$ (or reversed) may occure. Such equilibrium transitions are acompanied by the reversible jumps in the magnetization (see the broken vertical lines in Fig. 12($`b`$)). Thus, the equilibrium magnetization curve should be reversible and have a characteristic saw-like shape. Note also, that in the case of equilibrium transitions the magnetization should be positive for all external fields ($`4\pi M_\lambda >0`$), this corresponds to the diamagnetic susceptibility. \[The Meissner state with $`m=0`$ is always diamagnetic.\]

However, if the metastable states with $`m>0`$ are also admissible (i.e. those parts of curves in Fig. 12($`a`$), which lay to the left of the intersection points), then the magnetization may turn negative within some field interval ($`4\pi M_\lambda <0`$ for $`m1`$), this corresponds to the paramagnetic susceptibility. The magnetization curve in this case might be irreversible and display the hysteresis behavior.

Note also, that some of the magnetization curves in Fig.12($`b`$) ($`R_\lambda =2`$) have smooth ”tails” (for instance, at $`m=0`$), what indicates to the presence of the edge-suppression effect. For larger radii ($`R_\lambda >R_\mathrm{\Delta }4`$), instead of the inflection points, the reversible jumps in the magnetization may appear, which are analogous to those in Fig. 6($`a`$). Probably, such qualitative details may be seen in the experiments with mesoscopic samples.

As is evident from Fig. 12($`c`$), each of the superconducting $`m`$-states can exist only within a limited interval of fields, where $`\psi _{max}>0`$. When the field $`h_\xi `$ is increased, only the last state (with $`m=11`$) survives. In this state the superconducting region with $`\psi 0`$ is maximally pushed toward the cylinder surface, but the interior is almost normal (see Fig. 13 below). Evidently, the state with $`m=m_{max}`$ is just the state of the surface superconductivity . In a massive cylinder ($`R_\lambda 1`$) the last traces of the superconductivity may persist up to the field $`h_\xi =1.69`$ \[or, the third critical field $`H=H_{c3}=1.69H_{c2}`$\]. At smaller $`R_\lambda `$ the superconductivity may persist in the fields, which exceed $`H_{c3}`$ (see Fig. 1($`a`$)).

A peculiar field-compression effect is also present in Fig. 12. The first superconducting solution, which appears in a field-cooled regime, has maximal possible $`m`$ (the state of the surface-enhanced superconductivity). When the external field diminishes, this state becomes impossible, but new $`m1`$-state appears, in which the superconducting region is shifted toward the center, and simultaneously the edge-suppressed region begins to form. The state $`m=4`$, for instance, has both the region of edge-suppressed superconductivity (where the external field is practically not screened), and the region at some distance from the center, where the field of the giant vortex ($`m=4`$) is compressed and the superconductivity is also almost crushed (because at the center the order parameter is always small, $`\psi (\rho )\rho ^m1`$). (A flux-compression effect was discussed also in ).

As was mentioned, Fig. 12($`b`$) indicates to the possible paramagnetic susceptibility of the superconducting cylinder in the state with $`m>0`$. This topic also deserves a special comment.

The paramagnetic Meissner effect (or, Wohlleben effect) was initially observed in  and extensively discussed in the literature (see, for example, Refs. \[6,8,11,15,41-43\] and references therein). A number of controversial explanations of this effect were proposed, but there is as yet no clear understanding of its nature. This effect is observed, in particular, in a finite-dimension samples, in the field-cooled regime, it shows the signs of metastability, of the reentrant behavior and hysteresis . It is interesting, that all these features follow naturally from Eqs. (1)–(6) of the Ginzburg–Landau theory.

Consider, for instance, the state $`m=4`$ (which is distinguished, for clarity, by thick curves in Fig. 12). This state can not exist in small external fields $`h_\xi `$, because the vortex magnetic field ($`m=4,R_\lambda =2`$) is too strong for a mesoscopic sample, and can not be confined inside. However, the imposition of a finite external field prevents the field dissipation, the situation stabilizes and the state $`m=4`$ becomes possible again in a larger field (the reentrance effect). Evidently, it is more easy to occupy the state $`m=4`$, going down from the stronger fields $`h_\xi `$ (in a field-cooled regime). The negative part of the magnetization curve $`m=4`$ (the paramagnetic branch) corresponds to a metastable state of higher free-energy (see Fig. 12. ($`a`$)), than the state $`m=3`$. \[Here the state $`m=4`$ is defined as metastable, on the ground, that its free energy is less than the free energy of the state $`m=3`$. Such metastability means, that the irreversible transition to the lower laying state $`m=3`$ is possible. However, the state $`m=4`$ continues to be stable relative small perturbations of its form.\]

Note, that in the positive part of the magnetization curve the superconducting currents flow, mainly, to screen out the external field. \[This leads to the diamagnetic response in Fig. 12($`b`$).\] In the negative part of the magnetization curve the role of the vortex field is more important, and the superconducting screening currents flow, mainly, in the opposite direction. \[This means the sign reversion of the magnetic moment, because it is proportional to the current, $`𝐌[𝐣_s𝐫]𝑑v`$. The presence of positive and negative contributions to the total magnetic moment is evident also from Fig. 11($`f`$), where the current density may be positive or negative). At the point $`M_\lambda =0`$ the currents, which screen the internal and external fields, counterbalance each other. Thus, according to the Ginzburg–Landau theory, the paramagnetism present in Fig. 12, has purely electrodynamic explanation.

Now we shall demonstrate, how the order parameter vanishes in approaching the critical curves $`h_c^{(m)}>0`$ in Figs. 1 and 2 \[the curves $`h_c`$ are supplied with an additional index $`(m)`$\]. The order parameter is small in the vicinity of $`h_c^{(m)}`$. In Fig. 13 the degenerated solutions (i.e. those with $`\psi _{max}1`$) are shown for $`R_\lambda =2`$, $`\text{ϰ}=2`$ and different $`m`$. \[The curves in Fig. 13 are normalized to unity, with true values $`\psi _{max}1`$.\] In the origin of coordinates the order parameter behaves as $`\psi (\rho )\rho ^m`$, in agreement with the theory (see, for instance, ). At $`\rho 1`$ the order parameter is also small (here the edge-suppressed state $`e_m`$ is formed). For all degenerated solutions, depicted in Fig. 13 (with $`\psi _{max}1`$), the field $`B(r)H`$. In this case Eqs. (1),(2) can be linearized and reduced to one linear equation for the function $`\psi `$ . The solutions of this equation may be found analitically, in terms of the Kummer functions (see, for instance, ). \[Notice, however, that the normalization factor in the linear equation for $`\psi `$ remains arbitrary, and to find the true value $`\psi _{max}`$ it is necessary to solve the full system of nonlinear equations (1)–(6).\] The critical fields $`h_c^{(m)}>0`$ in Fig. 1, which we found by solving the full system of equations (1)–(6), may be obtained more easily on the base of a linear equation (when $`\psi _{max}1`$), as was done in . \[However, when approaching the phase boundaries $`h_c^{(m)}`$ for negative $`h_\xi <0`$, or to find the critical fields $`h_s^{(m)}`$, when the edge-suppressed state forms (see Figs. 1 and 2), but the order parameter $`\psi _{max}`$ remains finite, it is necessary to use the full system of equations (1)–(6).\]

It is evident from Fig. 13, that even the degenerated solutions (with $`\psi _{max}1`$) display the presence of the edge-suppresssion effect. If $`m=m_{max}`$, the order parameter is maximal at the cylinder surface. For $`m<m_{max}`$ the order parameter begins being suppressed at the boundary. When the field $`H`$ diminishes, only the states with relatively small $`m`$ survive \[with the inflexion points on the magnitization curves, see Figs. 12($`b,c`$)\], and finally the edge-suppressed solutions would form \[analogous to those, shown in Fig. 11, with the characteristic jumps, as in Figs. 3 and 6\].

Finally, it is of interest to note, that the transitions between states with different $`m`$ are also possible in type-I superconductors (see also ). This is demonstrated by Fig. 14. In Figs. 14($`a,b`$) the profiles of the self-consistent solutions $`\psi (r)`$ and $`b(r)`$ are shown for type-I superconductor ($`\text{ϰ}=0.5`$, $`R_\lambda =4`$) in the fields: $`h_\xi =1.4892`$ ($`m=0`$); $`h_\xi =2.0344`$ ($`m=1`$); $`h_\xi =2.1084`$ ($`m=2`$). These $`h_\xi `$-values are in the immediate vicinity of the corresponding critical curves $`h_c^{(m)}`$ in Fig. 1($`b`$): if $`h_\xi `$ is increased by $`110^4`$ \[i.e., if the critical curves in Fig. 1($`b`$) are crossed along the line $`R_\lambda =4`$\], the transitions to the normal state occures ($`\psi (r)=0`$, $`B(r)=H`$).

As is clear from Fig. 14($`a`$), in type-I superconductors the transition to the normal state is of the first order (by a jump from a finite value $`\psi _{max}^{}`$). However, the jump amplitude (i.e. the value $`\psi _{max}^{}`$ at the transition point) depends on parameters $`m`$, ϰ and $`R_\lambda `$ (or $`R_\xi =R_\lambda /\text{ϰ}`$). The jump amplitude, $`\psi _{max}^{}(R_\lambda )`$, is depicted in Fig. 14($`c`$) for $`\text{ϰ}=0.1`$ and $`\text{ϰ}=1/\sqrt{2}`$ in the case $`m=0`$. As follows from Fig. 14($`c`$), the value $`\psi _{max}^{}`$ diminishes, when $`R_\lambda `$ diminishes, with the jump amplitude vanishing at the point $`R_\mathrm{\Delta }`$. Thus, the points $`R_\mathrm{\Delta }`$ on the critical curves in Figs. 1($`b`$), 10($`a`$) and 14($`c`$) divide two regions of the phase transitions in the case of type-I superconductors: for larger $`R_\lambda >R_\mathrm{\Delta }`$ the first order phase transition to the normal state occures; for smaller $`R_\lambda <R_\mathrm{\Delta }`$ one has the second order phase transition. At $`R_\lambda =R_\mathrm{\Delta }`$ both transitions became indistinguishable, because the jump amplitude $`\psi _{max}^{}`$ vanish. \[Remind, that in type-II superconductors, with $`\text{ϰ}>1/\sqrt{2}`$, the phase transition to the normal state is always second order, regardless of the specimen size.\] The critical size $`R_\mathrm{\Delta }`$, which separates the regions of the first and second order phase transitions in the case of type-I superconductors, may also be found from different considerations . Fig. 14($`c`$) shows, that this size depends slightly on ϰ.

## IV Comments and conclusions

In conclusion, the main results of the present investigation are formulated.

The self-consistent solutions of a nonlinear system of Ginzburg–Landau equations, which correspond to the fixed values of the magnetic quantum number $`m`$ and have a cylinder symmetry (i.e. depend only on the radial coordinat $`r`$, but not on the azymuthal angle $`\phi `$), are studied.

It is found, that for type-II superconducting cylinder of finite radius ($`R>R_\mathrm{\Delta }`$), placed in sufficiently weak magnetic field $`H`$, the order parameter $`\psi (r)\mathrm{const}`$ near the cylinder surface, and behaves as $`\psi (r)r^m`$ near the center (the giant-vortex state). If the field is increased ($`H>H_s^{(m)}`$), the giant-vortex state changes abruptly and acquires new form. (A first order phase transition to the edge-suppressed state occures.) In this new state the order parameter is strongly suppressed near the cylinder edge, with the magnetic field penetrating practically without screening into some layer in the vicinity of the surface. At the center the behavior $`\psi (r)r^m`$ remains. This transition to the edge-suppressed form is acompanied by jumps in the magnetization and in the total magnetic flux, confined in the cylinder. The properties of such edge-suppressed solutions (which depend on the parameters $`m,H,\text{ϰ},T`$) are studied. It is shown, that the edge-suppressed states exist only in type-II superconductors, but the giant-vortex states are possible both in type-I and type-II superconductors. The phase boundaries are found, inside which the giant $`m`$-vortex solutions can exist. In particular, the characteristic radius of the cylinder is found, $`R_\mathrm{\Delta }`$, which separates the regions of the first and second order phase transitions between the superconducting and normal states.

The behavior of the magnetization ($`M`$) and of the total magnetic flux ($`\mathrm{\Phi }_1`$) as functions of the external field ($`H`$) is found, with account of the jumps, which acompany the formation of the edge-suppressed states. It is demonstrated, that some maximal value $`m=m_{max}`$ exist, for which the edge-suppressed state degenerates into the usual state of the surface superconductivity (with nucleation field $`HH_{c3}`$ ). The Gibbs free energies of different axially symmetrical superconducting states are compared, and the problem of paramagnetic effect, observed in mesoscopic samples, is shortly discussed. It is found, that the paramagnetic effect is possible in metastable vortex states and it may be attributed to the disbalance of the superconducting currents, which screen the external field $`H`$ and the field of the vortex from entering the sample interior.

As was pointed out earlier , according to the Ginzburg–Landau theory, two competing mechanisms of the external field penetration into the superconducting cylinder exist – in a form of vortices \[2–15,30\], and in a form of the edge-suppressed layer. The question, which of these mechanisms is more favorable energetically, may be answered only after the full comparison is made of the Gibbs free energies for the arbitrary multy-vortex configurations in a finite radius cylinder. This difficult theoretical problem (as well as the generalisation to other geometries) is left outside the scope of the present investigation. Some of the edge-suppressed states, described above, are metastable, but they may be observable in a specially arranged experiments. We beleve it is expedient to draw attention to this possibility.

## V Acknowledgments

We are gratefull to V. L. Ginzburg for the interest to this work and valuable discussions.

Figures captions

Fig. 1.$`(a)`$ – The phase diagram of type-II superconductor ($`\text{ϰ}=2`$). The numerals at the curves – the values of the magnetic quantum number $`m`$. The value $`h_\xi =1`$ corresponds to the critical field $`H_{c2}=\varphi _0/(2\pi \xi ^2)`$. The value $`h_\xi =1.69`$ corresponds to the maximal field $`H_{c3}=1.69H_{c2}`$, in which the surface superconductivity is still possible in macroscopic specimens. ($`b`$) – The phase diagram of type-I superconductor ($`\text{ϰ}=0.5`$, $`m=0,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2}`$). The arrows and the letter $`R_\mathrm{\Delta }`$ at the curves denote the points, where the first and second order phase transitions to the normal state became idistinguishable (see the text).

Fig. 2. ($`a`$) – The case $`m=0`$, $`\text{ϰ}=2`$. The phase curve $`h_c`$ divides superconducting ($`s`$) and normal ($`n`$) states of the cylinder. When the representation point $`(R_\lambda ,h_\xi )`$ crosses the line $`h_s`$ (the dotted line), a new edge-suppressed solution forms, with suppressed order parameter $`\psi _1`$ (see Fig. 3($`f`$)). No edge-suppressed states exist for $`R_\lambda <R_w`$. ($`d`$) – Analogously, for the case $`m=1`$, $`\text{ϰ}=2`$. When $`R_\lambda <R_0`$, the field $`h_\xi >0`$ stimulates the superconducting state, so the reentry superconductivity is possible. ($`b`$) and ($`e`$) – Analogously, for different $`\text{ϰ}1`$. ($`c`$) and ($`f`$) – The phase curves of type-I superconductor (for various $`\text{ϰ}1/\sqrt{2}`$). The edge-suppressed states in type-I superconductors are abscent.

Fig. 3. The case $`m=0`$, $`\text{ϰ}=2`$, $`R_\lambda =5`$. ($`a`$) – The cylinder magnetic moment, $`M_\lambda =M/H_\lambda `$. ($`b`$) –The average magnetic field in the cylinder, $`\overline{b}`$ (or, the total magnetic flux in the system, $`\varphi _1=\frac{1}{2}\rho _1^2\overline{b}`$ (8)). ($`c`$) – The normalized difference of free energies, $`\mathrm{\Delta }g`$, (10). ($`d`$) – The magnetic field at the cylinder axis, $`b_0`$. ($`e`$) – The maximal value of the order parameter, $`\psi _{max}`$. ($`f`$) – The order parameter at the cylinder surface, $`\psi _1`$. The letter $`v_0`$ denotes the Meissner state (vortex-free, $`m=0`$). The letter $`e_0`$ denotes the edge-suppressed state (vortex-free, $`m=0`$). The transition from one branch of the solution to another takes place between points 2,3, or 5,6. (The solution profiles for the case $`R_\lambda =5`$, $`m=0`$, $`\text{ϰ}=2`$ are depicted in , Fig.2.)

Fig. 4. Analogous to Fig. 3, but for the case $`m=0`$, $`\text{ϰ}=2`$, $`R_\lambda =3`$.

Fig. 5. Analogous to Fig. 3, but for the case $`m=0`$, $`\text{ϰ}=2`$, $`R_\lambda =1`$.

Fig. 6. The same, as in Fig. 3, but for $`m=1`$, $`\text{ϰ}=2`$, $`R_\lambda =5`$. The letter $`v_1`$ denotes the Meissner-type vortex state (with a vortex $`m=1`$ at the cylinder axis). The letter $`e_1`$ denotes the edge-suppressed state with the vortex ($`m=1`$) at the axis. (The solutions profiles in the case $`m=1`$, $`\text{ϰ}=2`$, $`R_\lambda =5`$, are shown below in Fig. 11.)

Fig. 7. Analogous to Fig. 3, but for the case $`m=1`$, $`\text{ϰ}=2`$, $`R_\lambda =3`$.

Fig. 8. Analogous to Fig. 3, but for the case $`m=0`$, $`\text{ϰ}=2`$, $`R_\lambda =1`$.

Fig. 9. The same, as in Fig. 3, but for type-I superconductor ($`\text{ϰ}=0.5`$, $`R_\lambda =5`$, $`m=0`$ and $`m=1`$). There are no edge-suppressed states and no ”tails” in the magnetization for type-I superconductors.

Fig. 10. ($`a`$) – The phase diagram for $`m=0`$ and different ϰ. The arrow and letter $`R_\mathrm{\Delta }`$ devide the regions of the first and second order phase transitions to the normal state in the case of type-I superconductors. For $`m=0`$ the curves with different ϰ \[see numbers in Fig. 10($`a`$)\] practically coincide. ($`b`$) – The phase diagram for $`m=1`$ and different ϰ. For $`R_\xi <R_{\xi min}`$ the superconducting state with the vortex ($`m=1`$) at the axis is impossible.

Fig. 11. The coordinate dependence of solutions for $`m=1`$, $`R_\lambda =5`$, $`\text{ϰ}=2`$ and different fields (see the values in Fig. 11($`a`$). ($`a`$) – The order parameter, $`\psi (\rho )`$; $`\rho =r/\lambda `$. ($`b`$) – The derivative, $`d\psi /d\rho `$. ($`c`$) – The second derivative, $`d^2\psi /d\rho ^2`$. ($`d`$) – The field potential, $`U(\rho )`$. ($`e`$) – The normalized magnetic field, $`b(\rho )`$. ($`f`$) – The normalized current, $`j(\rho )`$.

Fig. 12. ($`a`$) – The free energy, $`\mathrm{\Delta }g`$; ($`b`$) – the magnetic moment, $`(4\pi M_\lambda )`$ and ($`c`$) – the maximal value of the order parameter, $`\psi _1`$, for the cylinder with $`R_\lambda =2`$, $`\text{ϰ}=2`$ and various $`m`$, as functions of $`h_\xi `$.

Fig. 13. The order parameter profiles in a case of degenerated states ($`\psi _{max}1`$) for $`R_\lambda =2`$, $`\text{ϰ}=2`$ and various $`m`$. (All the curves represent the self-consistent solutions of Eqs. (1)–(6), but they are normalized to unity.) With $`m`$ increasing, the value $`\psi _{max}`$ approaches the boundary $`r=R`$. The last of the degenerated states ($`m=11`$), which is maximally pushed toward the boundary, corresponds to the state of the surface superconductivity .

Fig. 14. ($`a`$) – The order parameter profiles in a case of type-I superconductor ($`\text{ϰ}=0.5`$, $`R/\lambda =4`$). The external fields (see the numerals at the curves) are in the immediate vicinity of the critical boundaries $`h_c^{(m)}`$ in Fig. 1($`b`$). If the field $`h_\xi `$ is increased by $`110^4`$, the jump transition to the normal state occures. ($`b`$) – The corresponding magnetic field profiles. ($`c`$) – The jump amplitude, $`\psi _{max}^{}`$, as function of $`R_\lambda `$ for $`\text{ϰ}=0.1`$ and $`\text{ϰ}=1/\sqrt{2}`$. \[The jump amplitude $`\psi _{max}^{}`$ is the value of $`\psi _{max}`$ in a field, which just preceeds the first order phase transition to the normal state. The value $`R_\mathrm{\Delta }`$, where the jump amplitude vanish (see Figs. 1($`b`$), 10($`a`$), 14($`c`$)) depends on ϰ.\]