format: text

text:
# The Density Matrix Renormalization Group Method and Large-Scale Nuclear Shell-Model Calculations

## I Introduction

The nuclear shell model  is one of the most powerful approaches for a microscopic description of nuclear properties. The low–energy structure of a given nucleus is described in this approach by assuming an inert doubly–magic core and then diagonalizing the effective hamiltonian within an active valence space consisting of at most a few major shells. Despite the enormous truncation inherent in this approach, the shell-model method as just described can still only be applied in very limited nuclear regimes, namely for those nuclei with a sufficiently small number of active nucleons or a relatively low degeneracy of the valence shells that are retained. The largest calculations that have been reported to date are for the binding energies of nuclei in the $`fp`$–shell through $`{}_{}{}^{64}Zn`$ .

For heavier nuclei or nuclei farther from closed shells, one is forced to make further truncations in order to reduce the number of shell-model configurations to a manageable size. The most promising approach now in use is to truncate on the basis of Monte Carlo sampling. In this way, it has recently proven possible to extend the shell model beyond the $`fp`$–shell to describe the transition from spherical to deformed nuclei in the Barium isotopes .

Another attractive possibility is provided by the Density Matrix Renormalization Group (DMRG), a method that was initially developed and applied in the framework of low–dimensional quantum lattice systems . For a recent review, see ref. . A simplified version of this method was applied  to the two-level pairing model, showing its convergence properties. Subsequently, the method was extended to finite Fermi systems and applied to a pairing problem of relevance to ultrasmall superconducting grains . This new methodology, which is referred to as the particle-hole (or p-h) DMRG, was recently applied to a first test problem of relevance to nuclear structure . The application involved identical nucleons moving in a single large-j shell under the influence of a pairing plus quadrupole interaction with an additional single-particle energy term that split the shell into doubly-degenerate levels. Comparing with the results of exact diagonalization, it was shown that the method leads to extremely accurate results for the ground state and for low–lying excited states without ever requiring the diagonalization of very large matrices. Furthermore, even when the problem was not amenable to exact solution, the method was seen to exhibit rapid exponential convergence. All of this suggests that the DMRG method may indeed be a practical means of carrying out reliable large–scale shell-model calculations.

Before proceeding to realistic applications of the method, it is important to clarify its various ingredients and assumptions. Furthermore, it is critical to optimize how we apply the methodology for the subsequent more realistic and complex applications which will follow. With that in mind, we report here the results of a more thorough investigation of the recent application of the DMRG method to the one-orbit pairing plus quadrupole model. Several of the questions that were hinted at in ref.  are now addressed, including the feasibility of applying the method to very much larger model spaces.

The outline of the paper is as follows. In Section II, we review the basic features of the DMRG method and provide a fairly comprehensive description of the theoretical ingredients required for its implementation. In Section III, we describe the one-shell model that we will be using for the first tests of the DMRG method in nuclear physics. In Section IV we present the results of these test calculations for a problem of 10 identical nucleon restricted to move in a single $`j=25/2`$ shell. This problem can be solved exactly using the Lanczos algorithm, providing us with a good testing ground of the methodology including several variants not discussed in ref. . In Section V, we discuss the application of the method to a problem involving 40 nucleons in a $`j=99/2`$ orbit, a problem well beyond the limits of exact diagonalization. In Section VI, we summarize the principal conclusions of this work and outline some future directions for investigation.

## II Overview of the p-h DMRG procedure

### A Qualitative description

The basic idea of the p-h DMRG method is to systematically take into account the physics of all single–particle levels. The p-h nature of the method enters through a separation of the active single-particle levels into two sets that are separated by the Fermi surface; those above the Fermi surface are called the particle levels and those below are called the hole levels. The procedure begins by first taking into account the most important levels, namely those that are nearest to the Fermi surface, and then gradually including the others in subsequent iterations. At each step of the procedure, a truncation is implemented both in the space of particle states and in the space of hole states, so as to optimally take into account the effect of the most important states for each of these two subspaces of the problem. The calculation is carried out as a function of the number of particle and hole states that are maintained after each iteration, with the assumption that these numbers are the same. This parameter, which we call $`p`$, is gradually increased and the results are plotted against it. Prior experience from other applications of the methodology suggests that the results converge exponentially with $`p`$. Thus, once we achieve changes with increasing $`p`$ that are acceptably small we terminate the calculation. Assuming that the convergence is rapid, the method permits us to achieve an accurate description of the low-lying states of the system without ever having to diagonalize the large matrices that would arise in the absence of the DMRG truncation strategy.

### B A bit more detail

It is useful to now put the above qualitative remarks about the DMRG method into a more mathematical context.

The general problem that we wish to solve is a standard shell model problem of $`N_\nu `$ neutrons and $`N_\pi `$ protons (total $`N=N_\nu +N_\pi `$) interacting via a one- plus two-body hamiltonian, with the neutrons and protons each restricted to an active set of single-particle orbitals. For simplicity of presentation, we will focus in this subsection on just one type of particle, as is the case for the test applications we will report. In a later subsection, we will note those modifications required for neutron-proton systems.

We begin by splitting the set of multiply-degenerate single-particle levels into an ordered set of doubly-degenerate levels. This is most naturally done by carrying out an axially-symmetric deformed Hartree Fock (HF) calculation for the system. Each of the resulting HF levels is then doubly-degenerate. Going from the set of spherical single-particle levels to the axially-symmetric set of doubly-degenerate levels is schematically illustrated in figure 1 for a problem involving 8 identical nucleons in the $`fp`$ shell.

Each of these doubly-degenerate single-particle levels admits precisely four states. One is the state with no particles, one is the state with one particle having positive $`m`$ projection, one is the state with one particle having negative $`m`$ projection, and the last is the state containing two particles, with both positive and negative $`m`$.

This set of doubly-degenerate single-particle levels naturally splits into two smaller sets, as likewise illustrated in figure 1. Those levels above the (non-interacting) Fermi surface are the particle levels and those below the Fermi surface are the hole levels. Reiterating what was stated earlier, we will systematically take into account all of the levels of the problem, by gradually moving away from the Fermi surface.

In this basis, the shell model hamiltonian takes the form

$`H`$ $`=`$ $`{\displaystyle \underset{\alpha m}{}}ϵ_{\alpha m}a_{\alpha m}^{}a_{\alpha m}`$ (1)
$`+`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{\alpha _1m_1\alpha _2m_2\alpha _3m_3\alpha _4m_4}{}}\alpha _1m_1,\alpha _2m_2|V|\alpha _3m_3,\alpha _4m_4`$ (3)
$`\times a_{\alpha _1m_1}^{}a_{\alpha _2m_2}^{}a_{\alpha _4m_4}a_{\alpha _3m_3},`$
with both a single-particle energy term and a two-body interaction.

At this point, let us assume that we have already treated some of the levels of the problem, perhaps the first $`n`$ particle levels above the Fermi energy (forming a particle block) and the first $`n`$ hole levels below the Fermi energy (forming a hole block). Let us also assume that we have precisely $`p`$ particle states in the particle block and $`p`$ hole states in the hole block, i.e. the same number for each.

Furthermore, we will assume that we have calculated and stored the matrix elements of all possible sub-operators of the hamiltonian within the $`p`$-dimensional particle and hole blocks, respectively. For a one- plus two-body hamiltonian, the operators of relevance are: $`a^{}`$, $`a^{}a`$, $`a^{}a^{}`$, $`a^{}a^{}a`$ and $`a^{}a^{}aa`$, as all others can be directly calculated from these using appropriate hermitian adjoint relations.

We now enlarge both the particle and hole blocks, by adding one more (doubly-degenerate) level to each. Clearly the number of states in the enlarged particle and hole blocks will be $`4p`$. These states can be written as $`|i,j=|i|j`$, where $`i=1,\mathrm{},p`$ and $`j=1,\mathrm{},4`$.

What we would like to do is to truncate both the particle and hole blocks from these $`4p`$ states to the optimum $`p`$ states, precisely the same number that we had before adding the new level.

How do we accomplish this?

The first step is to recalculate the matrix elements of all of the hamiltonian sub-operators in the enlarged particle and hole blocks. To show how this is done, we will focus on one specific set of matrix elements, those of the operator $`a_{\alpha _1m_1}^{}a_{\alpha _2m_2}`$. Its matrix elements in the enlarged block can be readily expressed as

$`i,j|a_{\alpha _1m_1}^{}a_{\alpha _2m_2}|k,l`$ $`=i|a_{\alpha _1m_1}^{}a_{\alpha _2m_2}|k\delta _{j,l}`$ (7)
$`+j|a_{\alpha _1m_1}^{}a_{\alpha _2m_2}|l\delta _{i,k}`$
$`+()^{n_k}i|a_{\alpha _1m_1}^{}|kj|a_{\alpha _2m_2}|l`$
$`()^{n_k}i|a_{\alpha _2m_2}|kj|a_{\alpha _1m_1}^{}|l,`$
where $`n_i`$ is the number of particles in state $`|i`$. Depending on which of the two subspaces the operator indices refer (either the levels $`1n`$ or the new level $`n+1`$), a different term in the sum applies. For example, if both $`\alpha _1,m_1`$ and $`\alpha _2,m_2`$ refer to single-particle states from the first $`n`$ levels, it is the first term that applies.

The key point is that all terms in (LABEL:enlarge) involve a product of a matrix element from the previous $`p\times p`$ space (levels $`1n`$) and a matrix element for the level that is being added. The matrix element from the previous space was already calculated and stored, whereas the one from the new level is very simple and can be written down analytically.

At this point we have calculated and stored all of the particle and hole matrix elements in their enlarged blocks. The next step is to construct the space of particle-hole states, by forming products of states in the particle and hole blocks, viz: $`|i_P,j_H=|i_P|j_H`$. This is usually referred to as the superblock.

We only build states in the superblock with the correct number of total particles $`N`$ for the problem of interest. That means that we only consider states with the same number of particles and holes. Furthermore, assuming that the hamiltonian is rotationally invariant, it suffices to consider only those states for which the total z-projection of the angular momentum is zero.

Within this product space, we then build the hamiltonian matrix. Using the simplified notation $`\gamma =\{\alpha m\}`$, this can be expressed as

$`i_P,j_H|`$ $`H`$ $`|k_P,l_H={\displaystyle }_\gamma ϵ_\gamma [i_P|a_\gamma ^{}a_\gamma |k_P\delta _{j_H,l_H}`$ (9)
$`+j_H|a_\gamma ^{}a_\gamma |l_H\delta _{i_P,k_P}]`$ (10)
$`+{\displaystyle \frac{1}{4}}{\displaystyle \underset{\gamma _i,\gamma _j,\gamma _k,\gamma _l}{}}\gamma _i,\gamma _j|V|\gamma _k,\gamma _l`$ (20)
$`[\delta _{i_P,k_P}j_H|a_{\gamma _i}^{}a_{\gamma _j}^{}a_{\gamma _l}a_{\gamma _k}|l_H`$
$`+2()^{n_{k_P}}i_P|a_{\gamma _i}^{}|k_Pj_H|a_{\gamma _j}^{}a_{\gamma _l}a_{\gamma _k}|l_H`$
$`+2()^{n_{k_P}}i_P|a_{\gamma _l}|k_Pj_H|a_{\gamma _i}^{}a_{\gamma _j}^{}a_{\gamma _k}|l_H`$
$`+i_P|a_{\gamma _i}^{}a_{\gamma _j}^{}|k_Pj_H|a_{\gamma _l}a_{\gamma _k}|l_H`$
$`4i_P|a_{\gamma _i}^{}a_{\gamma _l}|k_Pj_H|a_{\gamma _j}^{}a_{\gamma _k}|l_H`$
$`+i_P|a_{\gamma _l}a_{\gamma _k}|k_Pj_H|a_{\gamma _i}^{}a_{\gamma _j}^{}|l_H`$
$`+2()^{n_{k_P}}i_P|a_{\gamma _j}^{}a_{\gamma _k}a_{\gamma _l}|k_Pj_H|a_{\gamma _i}^{}|l_H`$
$`+2()^{n_{k_P}}i_P|a_{\gamma _i}^{}a_{\gamma _j}^{}a_{\gamma _l}|k_Pj_H|a_{\gamma _k}|l_H`$
$`\delta _{j_H,l_H}i_P|a_{\gamma _i}^{}a_{\gamma _j}^{}a_{\gamma _l}a_{\gamma _k}|k_H].`$
Note that every term in eq. (20) involves a product of a matrix element in the particle space and a matrix element of the conjugate operator in the hole space. All of these matrix elements, for both particles and holes, have already been calculated and stored in their respective enlarged blocks.

Next we diagonalize the superblock hamiltonian, viz:,

$$H|\mathrm{\Psi }_k=E_k|\mathrm{\Psi }_k,$$
(21)
with

$$|\mathrm{\Psi }_k=\underset{i_P,j_P=1,4p}{}\mathrm{\Psi }_{ij}^{(k)}|i_P|j_H.$$
(22)
Note that both sums go over the $`4p`$ states in the respective enlarged particle and hole blocks.

We now wish to construct the optimal approximation to the ground state of the system (or to some low-lying set of states) that is achieved when we only retain $`p`$ particle states and $`p`$ hole states. By optimal, we will mean the approximation that maximizes the overlap between the truncated state and the exact ground state of the superblock ($`|\mathrm{\Psi }_1`$ in eq. (22)). For now, we will ask for the optimal description of the ground state; later we will discuss how to generalize this to several states.

We will perform the optimized truncation in two stages, first asking what is the optimal approximation when we truncate the particle states and then asking what is the optimal approximation when we truncate the hole states. For simplicity of notation, we will let $`|\mathrm{\Psi }=|\mathrm{\Psi }_1`$ and $`\mathrm{\Psi }_{ij}=\mathrm{\Psi }_{ij}^{(1)}`$.

Let us focus first on the particle block. We wish to find a subset of particle states ($`p`$ in number), such that the truncated ground state built up from these states has the largest possible overlap with the exact superblock ground state. To accomplish this, we first define the reduced density matrix for particles,

$$\rho _{i_Pi_P^{}}^P=\underset{j_H=1,4p}{}\mathrm{\Psi }_{i_Pj_H}\mathrm{\Psi }_{i_P^{}j_H}^{}.$$
(23)
The reduced density matrix separates into blocks according to the number of particles and the $`m`$ value.

We then diagonalize this $`4p\times 4p`$ matrix,

$$\rho ^P|u_P^\alpha >=\omega _\alpha ^P|u_P^\alpha >.$$
(24)
A given eigenvalue $`\omega _\alpha ^P`$ represents the probability of finding the particle state $`|u_P^\alpha >`$ in the full ground state wave function $`|\mathrm{\Psi }>`$. Thus, the optimal truncation in the sense described above corresponds to retaining the $`p`$ eigenvectors $`|u_P^\alpha >`$ that have the largest eigenvalues.

We then do exactly the same for the hole block. Namely, we construct the reduced density matrix for holes,

$$\rho _{j_Hj_H^{}}^H=\underset{i_P=1,4p}{}\mathrm{\Psi }_{i_Pj_H}\mathrm{\Psi }_{i_Pj_H^{}}^{},$$
(25)
diagonalize it, and then keep only the $`p`$ hole states with the largest eigenvalues.

At this point, we have identified the $`p`$ particle states and $`p`$ hole states that best approximate the full ground state of the coupled superblock. The final step before proceeding to the next level is to transform all of the matrix elements that we calculated in the $`4p`$-dimensional particle and hole spaces to the optimal $`p`$-dimensional truncated spaces.

As noted earlier, we need not target our optimization to the ground state only. We can target it to any set of states we wish. If, for example, we wanted to optimally describe the lowest $`L`$ eigenstates of the system, we would build mixed density matrices

$`\rho _{i_Pi_P^{}}^P`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{k=1,L}{}}{\displaystyle \underset{j_H}{}}\mathrm{\Psi }_{i_Pj_H}^{(k)^{}}\mathrm{\Psi }_{i_P^{}j_H}^{(k)}`$ (26)
$`\rho _{j_Hj_H^{}}^H`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{k=1,L}{}}{\displaystyle \underset{i_P}{}}\mathrm{\Psi }_{i_Pj_H}^{(k)^{}}\mathrm{\Psi }_{i_Pj_H^{}}^{(k)}`$ (27)
and use them to choose the $`p`$ most important particle and hole basis states to retain.

Once this series of steps has been implemented, we simply return to the point where the discussion began and continue iteratively. Namely we add the next levels, precisely as we did above, first calculating and storing all matrix elements in the enlarged particle and hole blocks, then constructing and diagonalizing the hamiltonian matrix in the superblock, then constructing and diagonalizing the reduced density matrices for particles and holes (depending on the states we wish to target) and then truncating to the most important $`p`$ states for particles and holes, respectively, based on the eigenvalues of the associated density matrices.

Were the number of particle and hole levels the same, we would simply proceed as above to treat all levels. We would initialize the iteration process by treating the two levels nearest to the Fermi surface, whose matrix elements can be calculated very simply. We would then continue to add levels as defined above. Truncation would not be necessary until the number of particle/hole states exceeds $`p`$. Prior to that point, it is not necessary to construct the hamiltonian matrix in the superblock, since it is only needed for the purpose of truncation. Note that in the absence of truncation, the number of particle and hole states after treating $`n`$ levels is $`4^n`$. Thus, until $`4^n>p`$, truncation is not required.

When there are different numbers of particle levels and hole levels, some modification to the above algorithm is needed. Assume for example that there are more particle levels than hole levels, so that the system is less than half full. We can continue to add levels two at a time, one for particles and one for holes, until we have exhausted all of the hole levels. From that point on, however, there are no hole levels to add. Thus, in subsequent iterations, we only add particle levels and only carry out the optimized truncation for those states.

At each iteration, there are in fact $`4p`$ states for particles and holes in the enlarged blocks required to build the superblock hamiltonian. It is only after truncation that it is reduced to $`p`$. Assuming that there are more particle levels than hole levels, we have a choice of how to proceed after the hole levels have been exhausted. One possibility is not to carry out a truncation of the hole levels at its last iteration, keeping $`4p`$ hole states in subsequent iterations. The other is to carry out a truncation at the last hole iteration and then keep only $`p`$ hole states subsequently. The former is more accurate; the latter reduces (perhaps significantly) the storage needs for the calculation. Considering that storage issues are likely to be especially important for subsequent applications of the methodology to more complex systems, it is useful to assess the accuracy lost if we keep only $`p`$ hole states after the hole levels have been exhausted. This is one of the issues we will address in the test applications to follow.

### C Neutrons and Protons

In the presence of both neutrons and protons, the same formalism applies, but with a few minor practical modifications.

* In a system of neutrons and protons, there are four distinct blocks – neutron particle, proton particle, neutron hole and proton hole.
* In a given iteration, we have a choice of how many and which levels to include. One possibility is to add a particle level and a hole level of the same type, switching between proton levels and neutron levels from one iteration to the next. The other possibility is to add four levels in an iteration – one for each of the four blocks. Both options can readily be implemented with the formalism of the previous subsection and both permit us to eventually include all active levels. We expect that the first approach, in which we add a particle and a hole level in each iteration, should suffice for heavy nuclei where the active neutrons and protons occupy different major shells. For lighter nuclei where they occupy the same major shells, it is most likely preferable to add all four levels at the same time.
* The superblock basis in which the hamiltonian is diagonalized consists of products of states of all four types – neutron particle, neutron hole, proton particle and proton hole. Only those states with the correct numbers of neutrons and protons and with total angular projection zero need be considered.
* The reduced density matrix for a given block requires a contraction over the states of the other three blocks.
* We need to consider hamiltonians that include single-particle energy terms for neutrons and protons as well as the neutron-neutron, proton-proton and neutron-proton interactions. The matrix elements for all of these terms can be readily obtained from the same basic like-particle operators as before ($`a^{}`$, $`a^{}a`$, $`a^{}a^{}`$, $`a^{}a^{}a`$ and $`a^{}a^{}aa`$).

## III The model

The first application of the p-h DMRG, as described in the previous section, was to a problem of a very large number of particles interacting via a pairing force and constrained to a set of equally-spaced doubly-degenerate single-particle levels . This problem, which is of direct relevance to the physics of ultrasmall superconducting grains, was solved at half filling for up to 400 levels. Despite the enormity of the model space, especially for a very large number of levels, the DMRG method was able to very accurately reproduce the ground state of the system, which could be obtained exactly using a method pioneered by Richardson. These results suggest that the p-h DMRG method can very accurately describe the properties of fermion systems that are dominated by a single collective degree of freedom.

Nuclei, on the other hand, are characterized by several competing collective degrees of freedom, most importantly those associated with the pairing and quadrupole fields. Thus, before applying the DMRG method with confidence to nuclear systems, it is necessary to first demonstrate that it continues to work well in the presence of such competing collective effects.

With that in mind, we report here a series of test calculations for systems that admit pairing plus quadrupole correlations on the same footing. The simplest such system is one involving identical nucleons restricted to a single level with large degeneracy and interacting via a pairing plus quadruple force. Unfortunately, the ($`0^+`$) ground state for an even number of particles occupying a single level invariably has equal population of all $`m`$ substates and thus such a problem does not admit a Fermi surface. Since a Fermi surface is an important characteristic of real Fermi systems, including nuclei, it is important to permit one in the test calculations. Thus, in the calculations to follow we supplement the pairing plus quadrupole interaction of our test model with a single-particle energy term that splits the multiply-degenerate large-j oprbit into a series of doubly-degenerate levels. The hamiltonian we use takes the form

$$H=\chi QQgP^{}Pϵ\underset{m}{}|m|a_{jm}^{}a_{jm}.$$
(28)
The last term splits the levels of the single-j shell into a set of equally spaced levels of “oblate” character, with the largest $`|m|`$ value lowest. Because of the last term, the hamiltonian is not in general rotationally invariant and thus its eigenstates do not have conserved angular momentum. Rotational invariance is of course recovered for $`ϵ=0`$.

In the calculations we will describe, we use a scaled version of the quadrupole-quadrupole strength $`\chi `$, defined according to

$$\chi =\frac{jQj^2}{2j+1}\stackrel{~}{\chi },$$
(29)
where $`\stackrel{~}{\chi }`$ is the usual $`QQ`$ coupling strength. With this definition, $`\chi `$ and $`g`$ have the same dimensions, which we will subsequently ignore.

## IV Results for 10 particles in a j=25/2 orbit

### A General remarks

In this section, we report the results of calculations for a system of 10 identical particles in a $`j=25/2`$ orbit. In this case, there are 13 active doubly-degenerate orbits, 5 for holes and 8 for particles. The Fermi surface lies between the $`m=17/2`$ hole level and the $`m=15/2`$ particle. For this problem, the size of the space that would be required for exact diagonalization (assuming that we only consider $`m=0`$ states) is 109,583, well within the limits of standard diagonalization routines. Thus, for this set of calculations, we can readily compare the results of calculations based on the DMRG methodology with the corresponding exact results, and thereby obtain important insight into the usefulness of the method. We will first report calculations for the standard scenario, in which (i) we only target the ground state in the DMRG truncation strategy, and (ii) we include the full complement of $`4p`$ hole states in the hole block after all hole levels have been treated. We will then present results for scenarios that systematically modify these two assumptions.

### B Results for standard scenario

Here we present results in which we only target the ground state and we include all $`4p`$ hole states in the hole block after the hole levels have been exhausted.

In our earlier paper, we presented results for two cases, one with $`\chi =1,g=0,ϵ=0.1`$ and one with $`\chi =1,g=0.5,ϵ=0.1`$. In both cases, we achieved accuracy for the ground state energy of roughly $`1`$ part in $`10^6`$ with reasonable values of $`p`$.

Obviously the case with $`g=0`$ involves pure quadrupole correlations. On the basis of HFB calculations that we have now performed, we conclude that the $`g=0.5`$ case likewise is dominated by quadrupole correlations, with pairing correlations minimal. Since it is our desire here to assess the DMRG method in the typical nuclear scenario in which different correlation effects compete, we will show throughout this section results for a further-enhanced pairing interaction, namely for $`\chi =1,g=0.1,ϵ=0.1`$. In this case, the HFB calculations show a well developed superconducting solution. This can be readily seen by noting that for this problem the ratio of the average pairing gap to the mean spacing between HF levels is $`0.39/0.18=2.1`$.

Figure 2 shows results for the DMRG correlation energies, defined as the gain in energy relative to Hartree Fock approximation. The exact and HFB values are

$$E_{corr}^{Exact}=0.70633;E_{corr}^{HFB}=0.20641$$
(30)
By a $`p`$ value of $`200`$, we obtain more than 99% of the exact correlation energy, a dramatic improvement over HFB approximation which only gets 28%.

For $`p=200`$ the largest superblock matrix that had to be diagonalized had a dimension of $`2,886`$. For $`p=40`$, the maximum dimension was $`232`$; nevertheless, for this value of $`p`$ we still achieve over 95% of the full correlation energy.

Table 1 shows results for the excitation energies of the lowest states of this same system. Despite the fact that we only targeted the ground state in the density matrix phase of these calculations, the agreement is as good for the excitation energies of these low-lying states as it is for the ground state energy. By $`p=200`$, all three excitation energies are reproduced to 1% or better.

In Figure 3, we show the exact, DMRG and HFB results for the occupation numbers of the particle and hole levels associated with the ground state solution. The DMRG results, which are shown for $`p=200`$, are in excellent agreement with the exact results for all levels with appreciable occupation. In contrast, the HFB calculation does not obtain sufficient depletion of the Fermi sea.

From these results, we conclude that the p-h DMRG method is able to describe with extreme accuracy the low-lying properties of complex many-body systems with competing collective features.

### C Dependence on the number of states targeted

What is the importance of including excited states in the targeting procedure? We address this in Figure 4, where we compare the results for the ground state energy and the energies of low-lying excited states when different sets of states are targeted. The solid triangles refer to results when only the ground state is targeted ($`L=1`$). The open circles are the results obtained when the lowest four states are all targeted in the optimization procedure ($`L=4`$). We also include in the figure the exact results.

The results are not unexpected. The ground state is better described when it is the only state targeted in the optimization procedure. Furthermore, in such a scenario a reasonable description of low-lying excited states is achieved. When several states are targeted, the agreement gets worse, albeit only marginally, for the ground state, but becomes better for the excited states. What is most important, however, is that once $`p`$ becomes large enough, there is no discernable difference between the two sets of results, both of which are also in excellent agreement with the exact eigen-energies.

### D Dependence on how we treat the hole block after all hole levels have been exhausted

Next we address the effect of maintaining all $`4p`$ hole states following treatment of the last hole level. This was addressed by carrying out the DMRG calculations in two ways, namely keeping only $`p`$ hole states after the last hole iteration or keeping $`4p`$. For a value of $`p=200`$, we find for the ground state energy a value of $`16.7331`$ when only $`p`$ states are retained and a value of $`16.7333`$ when all $`4p`$ are kept. The gain in energy through the improved treatment of the hole states is just $`1`$ part in $`10^6`$. On the basis of these results, the subsequent calculations we describe will be based on the simplifying, but highly justified, assumption that we maintain just $`p`$ hole states in those iterations that only add particle levels.

## V Results for 40 particles in a j=99/2 orbit

The excellent quality of the results obtained in the test calculations of the previous section has encouraged us to treat a much larger system, one for which exact diagonalization is not possible. Following an optimization of the DMRG calculational methodology, we are now able to treat systems significantly larger than were feasible in ref. . Here we report the largest calculation we have so far carried out – for a system of 40 particles occupying a $`j=99/2`$ orbit. In this case, the exact calculation would involve a hamiltonian matrix of dimension $`3.84\times 10^{25}`$, obviously much too large to treat without dramatic truncation. The parameters of the hamiltonian of eq. (28) that were used in this calculation are: $`\chi =1,g=0.1`$ and $`ϵ=0.2`$. The calculations targeted only the ground state and retained only $`p`$ hole states after treating the $`20`$ hole levels. \[The latter is especially important when treating extremely large problems such as this one.\] The results for the ground state correlation energy are shown in figure 7. By $`p=100`$, we obtain a ground state correlation energy of $`\mathrm{\hspace{0.17em}2.8994}`$. An exponential fit to the DMRG results, likewise indicated in the figure, give an asymptotic result for the ground state energy of $`\mathrm{\hspace{0.17em}2.8990}\pm 0.0003`$. This is close enough to our DMRG result with $`p=100`$ to suggest that we have achieved accuracy for the ground state correlation energy to 1 part in $`10^4`$.

We have also carried out an HFB calculation for this system, which confirms that the ground state is superconducting. The correlation energy achieved in this calculation is $`\mathrm{\hspace{0.17em}1.7902}`$, less than 62% of that estimated from the exponential fit to the DMRG results.

## VI Conclusion

In this paper, we describe the recently-developed particle-hole Density Matrix Renormalization Group method and report its test application to a problem involving many identical nucleons constrained to a single-j shell and subject to a hamiltonian with a pairing plus quadrupole interaction and a single-particle term that splits the shell into doublets. Earlier work along these lines was reported in ref. . Here we develop the formalism in greater detail and report calculations that address some important aspects of the method that needed to be clarified before the method could be meaningfully applied to more realistic nuclear systems.

In our view, the current calculations confirm and expand on the conclusions reported in ref.  regarding the use of the method in large-scale shell-model calculations. As long as there is a well defined Fermi surface in the problem, the method is able to produce extremely accurate results for the ground state energy of the system and for the energies of low–lying excited states as well, even in the presence of competing collective correlation effects. It is also able to accurately reproduce other (number-conserving) properties of the system, including for example occupation numbers. Based on an optimized computational methodology, we are now able to treat much larger problems than in ref. , and this bodes well for subsequent more-realistic applications of the methodology. In all cases we have studied, even those very large-scale problems we are now able to handle only because of the optimized methodology, we seem to be able to obtain accurate results while diagonalizing matrices of moderate dimensions. Critical to the success of the method is the rapid exponential convergence that it typically exhibits as a function of the number of particle and hole states maintained in a given calculation.

Clearly the next step is for us to further develop the method for use in realistic nuclear systems. This will mean including both neutrons and protons, several non-generate single-particle levels, and general nuclear hamiltonians. The framework for these extensions is spelled out in Section II and is currently in the process of being implemented.

###### Acknowledgements.

This work was supported in part by the National Science Foundation under grant # PHY-9970749, by the Spanish DGI under grant BFM2000-1320-C02-02, by NATO under grant PST.CLG.977000, by the Bulgarian Science Foundation under contract $`\mathrm{\Phi }905`$ and by the Bulgarian–Spanish Exchange Program under grant # 2001BG0009. One of the authors (SDD) would also like to acknowledge the partial support of a Fulbright Visiting Scholar Grant


text:
# 1 Introduction

## 1 Introduction

Buoyancy effects in unstably stratified mixing layers express themselves through gravity currents of heavy fluid which propagate in an ambient lighter fluid. These currents are encountered in numerous geophysical flows, industrial safety and environmental protection issues \[304haertel\] . During transition to turbulence a strong distortion of the separating interface between regions containing ‘heavy’ or ‘light’ fluid arises. The complexity of this interface will be used to monitor the progress of the mixing. We concentrate on the enhancement of surface-area and ‘surface-wrinkling’ of the separating interface as a result of gravity-effects.

Mixing in a turbulent flow may be visualized by a number of snapshots, e.g., in an animation. However, for a quantitative analysis, more specific measures are needed. Here we represent the separating interfaces as level-sets of a scalar field $`c`$ ($`0c1`$). We establish a close correspondence between the growth of surface-area and wrinkling, in particular at large buoyancy effects. We also show that this process can be simulated quite accurately using large-eddy simulation with dynamic subgrid modeling. However, the subgrid resolution, defined as the ratio between filter-width $`\mathrm{\Delta }`$ and grid-spacing $`h`$, should be sufficiently high to avoid contamination due to spatial discretization error effects.

## 2 Surface-area and wrinkling of level-sets

To specify the intricate structure of the separating interface we consider global properties of the level-set defined by $`S(t)=\{𝐱|c(𝐱,t)=1/2\}`$. The mixing-efficiency $`\eta `$ is defined in terms of the surface-area $`A`$ as:

$$\eta (t)=\frac{A(t)}{A(0)};A(t)=_{S(t)}𝑑S=_V𝑑𝐱\delta (c(𝐱,t)\frac{1}{2})|c(𝐱,t)|$$
(1)
Instead of integrating directly over the complex and possibly fragmented interface $`S(t)`$, we evaluate the volume integral in 1 , in which $`SV`$, using a locally exact quadrature method developed in  \[304jot2001\] . A further characterization of the geometry is obtained by considering the average curvature $`C`$ and wrinkling $`W`$:

$$C(t)=_{S(t)}𝑑S\kappa (𝐱,t);W(t)=_{S(t)}𝑑S|\kappa (𝐱,t)|$$
(2)
where the local curvature $`\kappa `$ is defined as

$$\kappa (𝐱,t)=𝐧;𝐧(𝐱)=\frac{c(𝐱)}{|c(𝐱)|}$$
(3)
The effect of gravity is represented by source terms in the momentum and energy equations \[304jot2001\] . Within the Boussinesq approximation we have:

$`_t\rho +_j(\rho u_j)=0`$
$`_t(\rho u_i)+_j(\rho u_iu_j)+_ip_j\sigma _{ij}=\left({\displaystyle \frac{\alpha \rho c}{Fr^2}}\right)\delta _{i2}`$
$`_te+_j((e+p)u_j)_j(\sigma _{ij}u_i)+_jq_j=\left({\displaystyle \frac{\alpha \rho c}{Fr^2}}\right)u_2`$ (4)
where $`\rho `$ denotes the ‘light’ fluid density, $`u_i`$ is the velocity in the $`x_i`$ direction, $`p`$ the pressure, $`e`$ the total energy density, $`q_j`$ the heat-flux and $`\sigma _{ij}=S_{ij}/Re`$ with $`S_{ij}`$ the rate of strain tensor and $`Re`$ the Reynolds number. The source terms on the right hand side represent gravity effects in the negative vertical (i.e. $`x_2`$) direction in which $`\alpha =(\rho _h\rho )/\rho `$ with $`\rho _h`$ the density of the ‘heavy’ fluid. The Froude number $`Fr=U^{}/\sqrt{g^{}L^{}}`$ with $`g^{}`$ the gravitational constant and $`U^{}`$, $`L^{}`$ a reference velocity and length scale respectively. We consider $`\alpha =1`$. Finally, $`c`$ denotes a scalar field which relates to the local density $`\stackrel{~}{\rho }`$ as $`\stackrel{~}{\rho }=(1c)\rho +c\rho _h`$. The scalar field $`c`$ is assumed to evolve according to

$$_t(\rho c)+_j(u_j\rho c)\frac{1}{Sc}_{jj}c=0$$
(5)
where $`Sc`$ denotes the Schmidt number. The scalar field $`c`$ is convected by the flow field $`u_j`$ but in turn affects the flow through the source terms in 4 . The initial state for $`c`$ is a step from a value 0 in the lower half to 1 in the upper half of the domain, representing a heavy fluid layer on top of a lighter one.

## 3 Buoyancy effects in turbulent mixing

To quantify buoyancy effects in turbulent mixing we consider the temporal mixing layer as documented in \[304vremanjfm\] . The flow starts from a perturbed laminar profile. A fourth order accurate finite volume discretization was used. Rapid transition to turbulence emerges, displaying helical pairings. At large buoyancy effects the separating interface becomes very distorted and protruding heavy/light fluid ‘fingers’ with a complex geometrical structure emerge, as shown in figure 1 (a).

The geometrical properties $`A`$, $`C`$ and $`W`$ of the separating interface are compared in figure 1 (b) for two different Froude numbers. The surface-area and wrinkling display approximately the same behavior, while the average surface-curvature remains quite limited, indicating that regions with negative and positive local curvature are about equally important in this flow.

In figure 2 we compare $`\eta `$ and $`W`$ for various $`Sc`$. An increase in $`Sc`$ reduces dissipation in 5 and results in enhanced mixing. After a laminar development up to $`t10`$, a nearly algebraic growth $`\eta t^\beta `$ emerges with $`\beta =\beta (Sc)`$. The wrinkling also develops algebraically. Initially this development is independent of $`Sc`$. Beyond $`t10`$ the algebraic growth-rate increases with $`Sc`$.

We also considered LES of this mixing flow. In figure 3 we compare predictions from DNS and LES with the dynamic mixed model. The filter-width $`\mathrm{\Delta }=L/16`$ with $`L`$ the side of the flow domain. The mixing-efficiency is slightly under-predicted on $`32^3`$ cells. The discrepancy with DNS is largely of numerical nature. By increasing the subgrid-resolution $`r=\mathrm{\Delta }/h`$ from 2 on $`32^3`$ to 4 on $`64^3`$ a strong improvement arises even though $`\mathrm{\Delta }`$ was kept fixed. The remaining errors are more representative of subgrid modeling deficiencies  \[304geurtsfroehlich2002\] .

## 4 Concluding remarks

Buoyancy can strongly increase turbulent mixing. A reasonable agreement between surface-area and wrinkling was observed. This may be relevant, e.g., in combustion modeling to capture flame-speed and intensity. Prediction of mixing-efficiency using LES with dynamic models, showed a large influence of spatial discretization errors. A subgrid-resolution $`r=4`$ is suggested.

\[ this page is intentionally left blank \]


text:
# 1 Introduction

## 1 Introduction

In superstring theories, there generally exist various dilaton and modulus fields. These fields (we call them ”moduli” ) are expected to acquire masses of the order of the gravitino mass $`m_{3/2}`$ through some non-perturbative effects of supersymmetry (SUSY) breaking. It is well known that the moduli cause serious cosmological problem  because they have only gravitationally suppressed interactions with other particles and hence have long lifetimes. The moduli with mass $`O(100)`$ GeV decay at the Big Bang Nucleosynthesis (BBN) epoch and spoil the success of the BBN by destroying the synthesized light elements, while moduli with lighter mass ($`_{}^<1`$ GeV) may overclose the universe or emit X($`\gamma `$)-rays giving too many contributions to the cosmic background radiation .

The mass of moduli ($`=`$ gravitino mass $`m_{3/2}`$) depends on models of SUSY breaking. In hidden sector models  the SUSY breaking in the hidden sector is mediated by gravitation and SUSY particles (squarks, sleptons, etc ) in the observable sector as well as gravitino obtain the mass of weak scale $`O(100)`$ GeV. On the other hand, in gauge-mediated SUSY breaking models , the gauge interactions mediate the SUSY breaking effects. In gauge-mediated SUSY breaking, the gravitino cannot acquires mass by the gauge interactions but only through gravitation. Thus, the mass of gravitino is much lighter ($`_{}^<1`$ GeV) than that in the hidden sector models. Here we consider the moduli problem in the gauge-mediated SUSY breaking models.

In order to avoid the moduli problem, we need some huge entropy production process by which the moduli density is diluted. So far the most successful mechanism for entropy production is “thermal inflation” proposed by Lyth and Stewart . The thermal inflation model in the gauge-mediated SUSY breaking models was intensively investigated in Refs.  where it was shown that the thermal inflation can solve the moduli problem. However, any process that dilutes the moduli also dilutes primordial baryon asymmetry of the universe. Thus, we must produce sufficiently large baryon asymmetry before the entropy production occurs. The promising candidate for mechanism of such efficient baryon number generation is the Affleck-Dine (AD) baryogenesis . In fact, it was shown in Ref.  that both the present baryon asymmetry and small moduli density can be explained by the thermal inflation and the Affleck-Dine mechanism for the gauge-mediated SUSY breaking models.

However, it has been found that the dynamics of the Affleck-Dine baryogenesis is complicated by the existence of Q balls . In the gauge-mediated SUSY breaking models the potential for the Affleck-Dine field becomes flat at large amplitudes. For such a flat potential the Q-ball formation is inevitable . In order to produce a large baryon number the initial amplitude of the AD field should be large, which also leads to the formation of Q balls with huge baryon($`=`$Q) number. Since large Q balls are stable, the baryon number may be confined in the form of Q balls and there may exist very small baryon asymmetry in the cosmic plasma, which means that the baryogenesis does not work.

Unstable Q balls can provide all the charges created before the charge trapping by the produced Q balls, but rather small amplitudes of the AD field are necessary for the Q balls to decay into the ordinary baryons, nucleons. Thus, sufficient baryon number is not created from the beginning.

Therefore we study the cosmological moduli problem and baryogenesis in the gauge-mediated SUSY breaking models taking into account the Q-ball formation. It is found that the Q balls seriously affect the Affleck-Dine baryogenesis and lower its efficiency***The detailed analysis is found in Ref. .

## 2 Moduli Problem

The modulus field $`\eta `$ obtain a mass of the order of the gravitino mass $`m_{3/2}`$. During the primordial inflation the modulus field is expected to sit at some minimum of the effective potential determined by the Kähler potential and the Hubble parameter. In general, the minimum during the inflation deviates from the true minimum of the moduli potential at low energies and the difference of the two minimum is considered to be of the order of the gravitational scale $`M(=2.4\times 10^{18}\text{GeV})`$. After the inflation, when the Hubble parameter becomes comparable to the mass of the modulus, the modulus field begins to roll down toward the true minimum and oscillates. Then, the modulus density is estimated as

$$\frac{\rho _{mod}}{s}\frac{1}{8}T_{RH}\left(\frac{\eta _0}{M}\right)^2,$$
(1)
where $`s`$ is the entropy density, $`T_{RH}`$ is the reheating temperature and $`\eta _0`$ is the initial amplitude of the modulus oscillation ($`\eta _0M`$). In deriving Eq.(1), we have assumed that the modulus mass is equal to $`m_{3/2}`$ and the reheating takes after the moduli field starts the oscillation. Since $`T_{RH}`$ should be higher than about 10 MeV to keep the success of the BBN, the moduli to entropy ratio is bounded from below,

$$\frac{\rho _{mod}}{s}_{}^>1.25\times 10^3\mathrm{GeV}.$$
(2)
The decay rate of the modulus is very small because it has only gravitationally suppressed interaction. The lifetime is roughly estimated as

$$\tau _\eta 10^{18}\mathrm{sec}\left(\frac{m_{3/2}}{100\mathrm{M}\mathrm{e}\mathrm{V}}\right)^3.$$
(3)
Thus, for $`m_{3/2}^{}{}_{}{}^{<}100`$ MeV, the lifetime is longer than the age of the universe and its present density much larger than the critical density which is given by $`\rho _c/s_0=3.6\times 10^9h^2\mathrm{GeV}`$, where $`h`$ is the present Hubble parameter in units of 100km/sec/Mpc and $`s_0(2.8\times 10^3`$ cm$`{}_{}{}^{3})`$ is the present entropy density. The modulus with larger mass (100 MeV $`_{}^<m_{3/2}^{}{}_{}{}^{<}`$ 1 GeV ) decays into photons whose flux exceeds the observed background X(or $`\gamma `$)-rays. Therefore the modulus is cosmological disaster and should be diluted by some entropy production process.

## 3 Affleck-Dine mechanism and Q-ball formation

In the Minimal Supersymmetric Standard Model (MSSM), there exist flat directions, along which there are no classical potentials. Since flat directions consist of squarks and/or sleptons, they carry baryon and/or lepton numbers, and can be identified as the Affleck-Dine field. These flat directions are lifted by SUSY breaking effects. In the gauge-mediated SUSY breaking models, the potential of a flat direction is parabolic at the origin, and almost flat beyond the messenger scale:

$$V_{gauge}\{\begin{array}{cc}m_\varphi ^2|\mathrm{\Phi }|^2\hfill & \mathrm{\Phi }M_S\hfill \\ M_F^4\mathrm{log}\frac{|\mathrm{\Phi }|^2}{M_S^2}\hfill & \mathrm{\Phi }M_S\hfill \end{array},$$
(4)
where $`M_S`$ is the messenger mass scale.

Since the gravity always exists, flat directions are also lifted by the gravity-mediated SUSY breaking effects:

$$V_{grav}=m_{3/2}^2\left[1+K\mathrm{log}\left(\frac{|\mathrm{\Phi }|^2}{M}\right)\right]|\mathrm{\Phi }|^2,$$
(5)
where $`K`$ is the numerical coefficient of the one-loop corrections. This term can be dominant only at high energy scales because of the small gravitino mass $`_{}^<O(1\text{GeV})`$.

In addition, there is a thermal effect on the potential, which appears at two-loop order, as pointed out in Ref. . This comes from the fact that the running of the gauge coupling $`g(T)`$ is modified by integrating out heavy particles which directly couples with the AD field. This contribution is given by

$$V_TT^4\mathrm{log}\frac{|\mathrm{\Phi }|^2}{T^2}.$$
(6)
The baryon number is usually created just after the AD field starts coherent rotation in the potential, and its number density $`n_B`$ is estimated as

$$n_B(t_{osc})\epsilon \omega \varphi _{osc}^2,$$
(7)
where $`\epsilon (_{}^<1)`$ is the ellipticity parameter, which represents the strongness of the A-term, and $`\omega `$ and $`\varphi _{osc}`$ are the angular velocity and the amplitude of the AD field at the beginning of the oscillation (rotation) in its effective potential.

Actually, however, the AD field feels spatial instabilities during its coherent oscillation, and deforms into nontopological solitons, Q balls . From numerical calculations , Q balls absorb almost all the baryon charges which the AD field obtains, and the typical charge is estimated as

$$Q\beta \left(\frac{\varphi _{osc}}{M_F}\right)^4,$$
(8)
where $`\beta 6\times 10^4`$. Consequently, the present baryon asymmetry should be explained by the charges which come out of the Q balls through the evaporation, diffusion, and decay of Q balls.

In the case of the unstable Q balls, they decay into nucleons and light scalar particles. Since the temperature at the BBN time is very low ($`1`$ MeV), Q balls cannot decay into light scalars. This rate is thus given by

$$\frac{dQ}{dt}_{}^<\frac{\omega ^3A}{192\pi ^2},$$
(9)
where $`A`$ is a surface area of the Q ball.

In the case of the stable Q balls, the evaporation is the only way to extract the baryon charges from Q balls. The total evaporated charge from the Q ball is estimated as ,

$$\mathrm{\Delta }Q10^{15}\left(\frac{m_\varphi }{\text{TeV}}\right)^{2/3}\left(\frac{M_F}{10^6\text{GeV}}\right)^{1/3}Q^{1/12}.$$
(10)
Hence the baryon number density is suppressed by the factor $`\mathrm{\Delta }Q/Q`$, in comparison with the case of no stable Q-ball production.

On the other hand, where $`V_{grav}`$ dominates the potential at larger scales, the gravity-mediation type Q balls (‘new’ type) are produced , if $`K`$ is negative, while, if $`K`$ is positive, it is not until the AD field enters $`V_{gauge}`$ dominant region that it feels instabilities, and the gauge-mediation type Q balls are produced (the delayed Q balls) . Notice that the sign of $`K`$ is in general indefinite in the gauge-mediated SUSY breaking models.

When the AD field starts to oscillate in the $`V_{grav}`$-dominant region, where $`H_{osc}\omega m_{3/2}`$, the baryon number is produced as $`n_B\epsilon m_{3/2}\varphi _{osc}`$. For the negative $`K`$ case, the ‘new-type’ Q balls are created, and its charge is written as

$$Q\stackrel{~}{\beta }\left(\frac{\varphi _{osc}}{m_{3/2}}\right)^2,$$
(11)
where $`\stackrel{~}{\beta }6\times 10^3`$. This type of the Q ball is also stable against the decay into nucleons, and the amount of the baryons in the present universe is explained by the charge evaporation from the Q balls. The charge evaporated from the Q ball is estimated as

$$\mathrm{\Delta }Q2.2\times 10^{20}\left(\frac{m_{3/2}}{100\text{keV}}\right)^{1/3}\left(\frac{m_\varphi }{\text{TeV}}\right)^{2/3}.$$
(12)
## 4 Baryogenesis and the moduli problem

As we mentioned in the Introduction, the late-time entropy production necessary for the dilution of the moduli also dilutes the baryon numbers created earlier very seriously, but the sufficient numbers could remain, if the Q-ball production is not taken into account. We will see that the Q-ball formation put very serious restriction on the efficiency of the AD baryogenesis, whether the produced Q balls are stable or not.

### 4.1 Stable Q balls

Since the baryon number is supplied only by the evaporation from stable Q balls, it should be suppressed by the factor $`\mathrm{\Delta }Q/Q`$, compared with no Q-ball formation. We will show that this fact considerably reduces the power of the AD baryogenesis.

#### 4.1.1 Gauge-mediation type Q balls when the zero-temperature potential is dominated

If the $`V_{gauge}`$ dominates over the thermal logarithmic potential $`V_T`$, $`M_{F}^{}{}_{}{}^{>}T_{osc}`$, so that the initial amplitude is constrained as $`\varphi _{osc}^{}{}_{}{}^{>}(T_{RH}/M_F)^2M`$. $`V_{gauge}`$ also dominates over $`V_{grav}`$, which leads to the condition $`\varphi _{osc}^{}{}_{}{}^{<}M_F^2/m_{3/2}`$. Combining these two equations, we have $`T_{RH}^{}{}_{}{}^{<}M_F^2/\sqrt{m_{3/2}M}`$.

From the stability ($`M_Q/Q<1`$ GeV) and survival condition ($`\mathrm{\Delta }Q<Q`$), the maximal baryon-to-entropy ratio is

$`Y_B`$ $`_{}^<`$ $`6.2\times 10^{27}\epsilon \left({\displaystyle \frac{\mathrm{\Omega }_{mod}h^2}{0.2}}\right)\left({\displaystyle \frac{m_\varphi }{\text{TeV}}}\right)^{2/3}.`$ (13)
Therefore, the baryon-to-entropy ratio is too small to explain the present value of the order $`10^{10}`$.

#### 4.1.2 Gauge-mediation type Q balls when the thermal logarithmic potential is dominated

In this case, $`V_T`$ is dominant over $`V_{gauge}`$, and the AD field starts to oscillate when $`H_{osc}T_{osc}^2/\varphi _{osc}`$. We should also use the charge of the formed Q ball as $`Q\beta (\varphi _{osc}/T_{osc})^4`$. Therefore, the fraction of the evaporated charge is written as

$`{\displaystyle \frac{\mathrm{\Delta }Q}{Q}}`$ $``$ $`8.9\times 10^{17}\left({\displaystyle \frac{m_\varphi }{\text{TeV}}}\right)^{2/3}\left({\displaystyle \frac{M_F}{10^6\text{GeV}}}\right)^{1/3}\left({\displaystyle \frac{T_{RH}}{M}}\right)^{11/3}\left({\displaystyle \frac{\varphi _{osc}}{M}}\right)^{11/2}.`$ (14)
Then the baryon-to-entropy ratio becomes

$`Y_B`$ $``$ $`5.3\times 10^{10}\epsilon \left({\displaystyle \frac{\mathrm{\Omega }_{mod}h^2}{0.2}}\right)\left({\displaystyle \frac{m_\varphi }{\text{TeV}}}\right)^{2/3}\left({\displaystyle \frac{M_F}{10^6\text{GeV}}}\right)^{1/3}\left({\displaystyle \frac{T_{RH}}{M}}\right)^{5/3}\left({\displaystyle \frac{\varphi _{osc}}{M}}\right)^{3/2}.`$ (15)
We plot the baryon-to-entropy ratio in the function of $`m_{3/2}`$ in Fig. 1. We can thus marginally explain the present value ($`Y_{B}^{}{}_{}{}^{>}10^{11}`$) for $`m_{3/2}^{}{}_{}{}^{<}100`$ keV, but the reheating temperature is $`2.0\times 10^{13}`$ GeV, so high for natural inflation models to provide. It is noticed that $`Y_B`$ is maximized at $`\mathrm{\Delta }QQ`$.

#### 4.1.3 Delayed Q balls when the zero-temperature potential is dominated

Both the AD field and moduli fields start to oscillate when $`H_{osc}m_{3/2}`$, the ratio of the baryon number and the energy density of the moduli stays constant to the present, the baryon number is given by $`n_B\epsilon m_{3/2}\varphi _{osc}^2\mathrm{\Delta }Q/Q`$, and the baryon-to-entropy ratio becomes

$`Y_B`$ $``$ $`2.8\times 10^{24}\epsilon \left({\displaystyle \frac{\mathrm{\Omega }_{mod}h^2}{0.2}}\right)\left({\displaystyle \frac{m_\varphi }{\text{TeV}}}\right)^{2/3}\left({\displaystyle \frac{m_{3/2}}{100\text{keV}}}\right)^{8/3}\left({\displaystyle \frac{M_F}{10^6\text{GeV}}}\right)^4\left({\displaystyle \frac{\varphi _{osc}}{M}}\right)^2.`$ (16)
Figure 2 shows the maximum value of the baryon-to-entropy ratio. As can be seen, this scenario is marginally successful ($`Y_B10^{11}`$) only for $`m_{3/2}200`$ keV. Notice that, the value of $`Y_B`$ is the same as that in the thermal logarithmic potential, to be considered in the next subsection, for $`m_{3/2}^{}{}_{}{}^{<}0.16`$ GeV.

#### 4.1.4 Delayed Q balls when the thermal logarithmic potential is dominated

In the case that $`V_T`$ is dominant over $`V_{gauge}`$, the maximum value of the baryon-to-entropy ratio is also shown in Fig. 2. We can see that $`Y_B`$ is marginally enough ($`Y_B10^{11}`$) only for $`m_{3/2}200`$ keV. Notice that the largest possible value of $`Y_B`$ is achieved when $`\mathrm{\Delta }QQ`$.

#### 4.1.5 New type Q balls

From Eqs.(11) and (12), we have

$$\frac{\mathrm{\Delta }Q}{Q}6.1\times 10^{23}\left(\frac{m_\varphi }{\text{TeV}}\right)^{2/3}\left(\frac{m_{3/2}}{100\text{keV}}\right)^{5/3}\left(\frac{\varphi _{osc}}{M}\right)^2.$$
(17)
This leads to the baryon-to-entropy ratio as

$`Y_B`$ $`_{}^<`$ $`8.8\times 10^{28}\epsilon \left({\displaystyle \frac{\mathrm{\Omega }_{mod}h^2}{0.2}}\right)\left({\displaystyle \frac{m_\varphi }{\text{TeV}}}\right)^{2/3}\left({\displaystyle \frac{m_{3/2}}{100\text{keV}}}\right)^{2/3},`$ (18)
which is too small to explain the present value.

### 4.2 Unstable Q balls

In the case of the unstable Q balls which decay into nucleons, it may destroy light elements synthesized at the BBN. Thus, a new constraint which the Q ball should decay before the BBN ($`1`$ sec), must be imposed.

#### 4.2.1 Gauge-mediation type Q balls when the zero-temperature potential is dominated

There are several condition to be imposed. The first is the condition that the Q ball is unstable, given by

$$\frac{\varphi _{osc}}{M}_{}^<2.6\times 10^6\left(\frac{M_F}{10^6\text{GeV}}\right)^2.$$
(19)
Second, the decay of the Q ball must be completed until the BBN, otherwise it would spoil the success of the BBN, so that the life time $`\tau _Q`$ should be less than about 1 sec. Thus, we have the following constraint:

$$\frac{\varphi _{osc}}{M}_{}^<1.0\times 10^6\left(\frac{M_F}{10^6\text{GeV}}\right)^{6/5}.$$
(20)
The largest possible value of $`Y_B`$ is obtained by the conditions Eq.(20) and $`M_{F}^{}{}_{}{}^{<}(m_{3/2}M)^{1/2}`$, and will be written as

$`Y_B`$ $`_{}^<`$ $`2.8\times 10^{19}\epsilon \left({\displaystyle \frac{\mathrm{\Omega }_{mod}h^2}{0.2}}\right)\left({\displaystyle \frac{m_{3/2}}{100\text{keV}}}\right)^{4/5}.`$ (21)
This is thus too small to explain the present value $`10^{10}`$.

#### 4.2.2 Gauge-mediation type Q balls when the thermal logarithmic potential is dominated

Since the Q-ball charge is expressed as $`Q\beta (\varphi _{osc}/M)^6(T_{RH}/M)^4`$, the unstable condition, $`M_FQ_{}^{1/4}{}_{}{}^{>}1`$ GeV, is given by

$$\frac{\varphi _{osc}}{M}_{}^<3.4\times 10^4\left(\frac{M_F}{10^6\text{GeV}}\right)^{2/3}\left(\frac{T_{RH}}{M}\right)^{2/3},$$
(22)
while the lifetime condition that the Q ball decays before the BBN time ($`1`$ sec), is written as

$$\frac{\varphi _{osc}}{M}_{}^<4.0\times 10^3\left(\frac{M_F}{10^6\text{GeV}}\right)^{2/15}\left(\frac{T_{RH}}{M}\right)^{2/3}.$$
(23)
From the above conditions, we have the baryon-to-entropy ratio as

$$Y_{B}^{}{}_{}{}^{<}8.7\times 10^{11}\epsilon \left(\frac{\mathrm{\Omega }_{mod}h^2}{0.2}\right)\left(\frac{m_{3/2}}{100\text{keV}}\right)^{2/5}.$$
(24)
We plot the largest possible value of $`Y_B`$ \[Eq.(24)\] in the function of $`m_{3/2}`$ in Fig. 3. As can be seen, we can explain the present value in the range $`m_{3/2}=10100`$ keV. However, the reheating temperature should be as high as $`5.6\times 10^{12}`$ GeV, which may be rather higher for the actual inflation models.

#### 4.2.3 Delayed Q balls

Firstly, it is found that the unstable condition can be expressed as $`m_{3/2}^{}{}_{}{}^{>}0.16`$ GeV both when $`V_{T}^{}{}_{}{}^{<}V_{gauge}`$ and when $`V_{T}^{}{}_{}{}^{>}V_{gauge}`$. Therefore, in either case, we can estimate the baryon-to-entropy ratio as

$$Y_{B}^{}{}_{}{}^{<}4.5\times 10^{17}\epsilon \left(\frac{\mathrm{\Omega }_{mod}h^2}{10^9}\right)\left(\frac{m_{3/2}}{0.16\text{GeV}}\right)^1\left(\frac{\varphi _{osc}}{M}\right)^2,$$
(25)
where $`\mathrm{\Omega }_{mod}h_{}^{2}{}_{}{}^{<}10^9`$ for $`m_{3/2}^{}{}_{}{}^{>}0.16`$ GeV from the X($`\gamma `$)-ray constraint. This is too small to explain the present value, even if $`\varphi _{osc}M`$.

## 5 Conclusion

We have investigate the possibility of the AD baryogenesis in the gauge-mediated SUSY breaking scenario, while evading the cosmological moduli problem by the late-time entropy production. In all the cases, the Q-ball formation makes the efficiency of the baryon number production considerably be diminish. In the zero-temperature potential $`V_{gauge}`$, whether the produced Q balls are stable or not, the largest possible baryon-to-entropy ratio is too small to explain the present value. This completely kills the successful situations considered in Ref. .

We have also found that there are some marginally successful situations when we take into account of the thermal effects on the effective potential of the AD field. However, these successful situations require very high reheating temperatures such as $`10^{12}10^{16}`$ GeV, which may be impossible to achieve in the actual inflation models.

In the delayed Q-ball formation case, we have found that enough baryon-to-entropy ratio can be created, even in the zero-temperature potential $`V_{gauge}`$ is dominant over the thermal logarithmic potential $`V_T`$. It might be the unique solution for the AD baryogenesis with solving the cosmological moduli problem, although the scale of $`M_F`$ is rather low.

In addition, successful situations above need the following conditions: the large initial amplitude such as $`\varphi _{osc}M`$, and $`\epsilon 1`$. This is realized if the A-terms, which make the AD field rotate in the effective potential, originate from some Kähler potential. Then, $`\epsilon (\varphi /M)^\gamma 1`$ for $`\varphi M`$, where $`\gamma >2`$.

ACKNOWLEDGEMENTS

This work is supported by the Grant-in-Aid for Scientific Research from the Ministry of Education, Science, Sports, and Culture of Japan, Priority Area “Supersymmetry and Unified Theory of Elementary Particles” (No. 707).


text:
# Vortex Interactions and Thermally Induced Crossover from Type-I to Type-II Superconductivity

## I Introduction

The nature of the phase-transition in systems of a scalar matter field coupled to a massless gauge-field has a long history in condensed matter physics, dating at least back to the introduction of the Ginzburg-Landau (GL) theory of superconductivity Ginzburg:1950 . At the mean-field level, ignoring spatial variations in gauge fields as well as matter fields leads to the prediction of a second order phase transition in the model, with classical mean-field exponents for all values of the GL parameter $`\kappa `$. The first attempt to seriously consider the role of fluctuations on the order of the metal-to-superconductor transition was made by Halperin, Lubensky, and Ma Halperin:1974 , who found that ignoring matter-field fluctuations entirely, and treating gauge-field fluctuations exactly, resulted in a permanent first order transition for all values of $`\kappa `$, since the gauge-field-fluctuations produced an extra term $`|\varphi |^3`$ in the matter field sector of the theory in three spatial dimensions, where the complex matter field is denoted by $`\varphi `$ and represents the condensate order parameter. (In the context of particle physics, Coleman and Weinberg Coleman:1973  studied the equivalent problem of spontaneous symmetry breaking due to radiative corrections in the Abelian Higgs model in four space-time dimensions, finding the additional term $`\varphi ^4\mathrm{ln}(\varphi ^2/\varphi _0^2)`$, where the real matter field is denoted $`\varphi `$ and represents a scalar meson.) Note\_order\_from\_disorder . Subsequently, Dasgupta and HalperinDasgupta:1981  found, using duality arguments in conjunction with Monte-Carlo simulations, that when gauge-field fluctuations and phase fluctuations of the scalar matter field are taken into account, but amplitude fluctuations are ignored, the phase transition is permanently second orderDasgupta:1981 . Bartholomew Bartholomew:1983  then reported results from Monte-Carlo simulations for the case when also amplitude fluctuations are taken into account, concluding that the phase transition changes from first to second order at a particular value of the GL parameter $`\kappa 0.4/\sqrt{2}`$. As far as this numerical value is concerned, note that the problem of finding a tricritical  value $`\kappa _{\mathrm{tri}}`$ separating first and second order transitions is extremely demanding even by present day supercomputing standards Mo:2001  (see below). Using ingenious duality arguments, Kleinert Kleinert:1982  obtained that the change from first to second order transition should occur at $`\kappa 0.8/\sqrt{2}`$. The value of $`\kappa `$ that separates a first order (discontinuous) transition from a second order (continuous) one, defines a tricritical point Note\_tricritical , and will hereafter be denoted $`\kappa _{\text{tri}}`$. Note that to obtain the above result, it is necessary to allow for amplitude fluctuations in the superconducting order parameter, which become important for small to intermediate values of $`\kappa `$, but are totally negligible in the extreme type-II regime $`\kappa >>1`$.

The critical properties of a superconductor may be investigated at the phenomenological level by the GL model of a complex scalar matter field $`\varphi `$ coupled to a fluctuating massless gauge field $`𝐀`$. It is this feature of the gauge field that makes the GL model so difficult to access by the standard techniques employing the renormalisation group Halperin:1974 ; Note . The GL model in $`d`$ spatial dimensions is defined by the functional integral

$`Z`$ $`={\displaystyle }𝒟A_\nu 𝒟\varphi \mathrm{exp}[{\displaystyle }d^dx[{\displaystyle \frac{1}{4}}F_{\mu \nu }^2+|(_\nu +iqA_\nu )\varphi |^2+m^2|\varphi |^2+\lambda |\varphi |^4]`$ (1)
where $`F_{\mu \nu }=_\mu A_\nu _\nu A_\mu `$, $`q`$ is the charge coupling the condensate matter field $`\varphi `$ to the fluctuating gauge field $`A_\mu `$, $`\lambda `$ is a self coupling, and $`m^2`$ is a mass parameter which changes sign at the mean field critical temperature. When all dimensionful quantities are expressed in powers of the scale represented by $`q^2`$, the GL model may be formulated in terms of the two dimensionless parameters $`y=m^2/q^4`$ and $`x=\lambda /q^2`$. In this case, $`y`$ is temperature like and drives the system through a phase transition, and $`x=\kappa ^2`$ is the square of the Ginzburg-Landau parameter. Depending on the value of $`x`$, the transition is either first order for $`x<x_{\text{tri}}`$, or continuous for $`x>x_{\text{tri}}`$ Kleinert:1982 ; Bartholomew:1983 ; Mo:2001 .

In a recent paperMo:2001 , we have determined $`x_{\text{tri}}=0.295\pm 0.025`$. This corresponds to a tricritical value of the GL parameter $`\kappa _{\text{tri}}=(0.76\pm 0.04)/\sqrt{2}`$, in rather remarkable agreement with the results of Ref. Kleinert:1982, . Moreover in Ref. Mo:2001,  it was also argued that this value of $`x`$ or $`\kappa `$ is the demarkation value which separates type-I and type-II superconductivity, rather than the classical mean-field value $`\kappa =1/\sqrt{2}`$. The connection can be made when one realizes that criticality at the metal-to-superconductor transition requires that topological defects of the matter field in the form of vortex loops are stable. On the other hand, there is a connection between critical exponents and geometric properties of a tangle of such vortex loopsHove:2000 . The fractal dimension $`D_H`$ of the vortex-loop tangle is connected to the anomalous scaling dimension $`\eta _\varphi `$ Nguyen:1999  of the matter field in a field theory of the vortex-loop gas, a theory dual to the original GL theory Tesanovic:1999 ; Nguyen:1999  by the relation $`D_H+\eta _\varphi =2`$. Since the anomalous scaling dimension is connected to the order parameter exponent $`\beta `$ of the dual matter field by the relation $`2\beta =\nu (d2+\eta _\varphi )`$ Nguyen:1999 ; Hove:2000 , it follows that a collapse of the vortex-loop tangle implies $`D_H=d`$ and hence $`\beta =0`$ indicative of a first order transition. Here, $`d`$ is the spatial dimension of the system. Now, a collapse of the tangle in turn implies an effective attraction between vortices, or type-I behavior. On the other hand, a stable vortex-loop tangle at the critical point, with fractal dimension $`D_H<d`$, implies first of all type-II behavior, but also $`\eta _\varphi >2d`$ and $`\beta >0`$, and hence a second order transition.

The above assertion, that the tricritical value of $`\kappa `$ separates first order from second order metal-to-superconductor transition, and moreover also separates type-I from type-II behavior when the system is on the phase-transition line $`y_c(x)`$, is in contrast to the conventional wisdom that type-I and type-II superconductivity is separated by $`x=0.5`$. Based on the above arguments, we have proposed the phase diagram shown in Fig. 1 of Ref. Mo:2001,  which contains a new line separating type-I and type-II superconductivity. The shape of this line was inferred from the observation that far from the phase transition, mean-field estimates of the boundary between type-I and type-II should be precise, and hence this boundary should asymptotically approach $`x=0.5`$ from below as the temperature is reduced.

It is the purpose of this paper to show directly, by computing the effective thermally renormalized interaction between vortices via large-scale Monte-Carlo simulations, that this quantity changes from being repulsive to attractive in the intermediate regime $`\kappa [0.76,1]/\sqrt{2}`$. Since the sign of the vortex-interaction is the microscopic diagnostics, in terms of vortex degrees of freedom, for distinguishing type-I from type-II superconductivity, the large-scale simulations we present in this paper confirm the above conjectures and plausibility arguments of Ref. Mo:2001, .

In an external field the GL model has classical solutions in terms of Abrikosov flux tubesAbrikosov:1957b , or Nielsen–Olesen vorticesOlesen:1973 , and the concept of type-I versus type-II superconductivity is based on the interaction between these vortices. For type-I superconductors they attract each other, whereas for type-II superconductors the interaction is repulsive. AbrikosovAbrikosov:1957b  showed that at the *mean field* level type-I and type-II superconductors are separated at $`\kappa =1/\sqrt{2}`$. We will refer to the value of $`\kappa `$ separating type-I from type-II behavior at $`\kappa _{\mathrm{I}/\mathrm{II}}`$, which we find varies with $`y`$. It is not a sharply defined quantity, since it represents a crossover line. The exception is at $`y=y_c,x=x_{\text{tri}}`$, where $`\kappa _{\mathrm{I}/\mathrm{II}}=\kappa _{\text{tri}}`$. Elaborate calculations of vortex interactions have been carried out Hartmann:1966 ; Kramer:1971 ; Bogomolnyi:1976 ; Bogomolnyi:1976b ; Jacobs:1979 , but none of these approaches take thermal fluctuations into account. A recent overview of superconductors with $`\kappa `$ close to $`1/\sqrt{2}`$ can be found in Ref. Lukyanchuk:2001, , see also Ref. KleinNog, .

Superconductors with $`\kappa 1/\sqrt{2}`$ were studied extensively in the 1960s and 1970sLukyanchuk:2001 , and in particular measurements on the metals Ta and Nb demonstrated that the notion of a *temperature independent* value of $`\kappa _{\mathrm{I}/\mathrm{II}}`$ was incorrectAuer:1973 . At the time, this was explained with a mean-field theory involving three GL parameters Maki:1969 . Thermal fluctuations, not addressed at the mean-field level, offer an alternative and above all simpler explanation for the observations of crossovers from type-I to type-II behavior in one and the same compound as the temperature is increased.

We have performed large scale Monte Carlo (MC) simulations on the lattice version of Eq. 1, with two vortices penetrating the sample in the $`\widehat{z}`$ direction. By measuring the interaction between these two vortices we have determined the value of $`\kappa _{\mathrm{I}/\mathrm{II}}`$, in particular how this value is affected by thermal fluctuations close to the critical point.

## II Model, simulations and results

To perform simulations on Eq. 1, we have defined a discrete version as followsKajantie:1999 :

$`Z`$ $`={\displaystyle 𝒟𝜶𝒟\psi \mathrm{exp}(S[𝜶,\psi ])}`$
$`S[𝜶,\psi ]`$ $`=\beta _G{\displaystyle \underset{𝐱,i<j}{}}{\displaystyle \frac{1}{2}}\alpha _{ij}(𝐱)^2{\displaystyle \frac{2}{\beta _G}}{\displaystyle \underset{𝐱,\widehat{ı}}{}}\mathrm{𝖱𝖾}\left[\psi ^{}(𝐱)e^{i\alpha _i(𝐱)}\psi (𝐱+\widehat{ı})\right]+\beta _2{\displaystyle \underset{𝐱}{}}\psi ^{}(𝐱)\psi (𝐱)+{\displaystyle \frac{x}{\beta _G^3}}{\displaystyle \underset{𝐱}{}}\left[\psi ^{}(𝐱)\psi (𝐱)\right]^2.`$ (2)
In Eq. II $`\alpha _i(𝐱)=aqA_i(𝐱)`$ and $`\alpha _{ij}=\alpha _i(𝐱)+\alpha _j(𝐱+\widehat{ı})\alpha _j(𝐱)\alpha _j(𝐱+\widehat{ȷ})`$. $`\beta _G`$ and $`\beta _2`$ are related to the continuum parameters $`x`$ and $`y`$ and the lattice constant $`a`$,

$`\beta _G`$ $`={\displaystyle \frac{1}{aq^2}}`$ (3)
$`\beta _2`$ $`={\displaystyle \frac{1}{\beta _G}}\left[6+{\displaystyle \frac{y}{\beta _G^2}}{\displaystyle \frac{3.1759115\left(1+2x\right)}{2\pi \beta _G}}{\displaystyle \frac{\left(4+8x8x^2\right)\left(\mathrm{ln}6\beta _G+0.09\right)1.1+4.6x}{16\pi ^2\beta _G^2}}\right].`$ (4)
Note that $`\beta _2`$ contains the effect of ultraviolet renormalization in the continuum limit when the lattice constant $`a0`$Laine:1998 ; Mo:2001 . The model Eq. II is defined on a numerical grid of size $`N_x\times N_y\times N_z`$, corresponding to a *physical* size of $`L_x\times L_y\times L_z`$, with $`L_i=N_ia`$. All our simulations have been on cubic systems with $`\beta _G=1`$.

To impose an external magnetic fieldKajantie:1999 , we modify the action Eq. II by changing the field energy along one stack of plaquettes located at $`x_0,y_0`$ in such a way that the action is minimized for $`\alpha _{12}(x_0,y_0,z)=2\pi n`$ instead of zero, corresponding to forcing a number of $`n`$ flux quanta through the system. Hence, the action $`S[\alpha ,\psi ;n]`$ for $`n`$ flux-quanta forced through the system, is given by

$$S[\alpha ,\psi ;n]=S[\alpha ,\psi ;0]+\underset{z}{}\left(2\pi n\alpha _{12}(x_0,y_0,z)+2\pi ^2n^2\right).$$
(5)
The second term in Eq. 5 corresponds to forcing a flux $`\mathrm{\Phi }_B`$ through the lattice in the negative $`z`$-direction

$$\frac{\alpha _{12}(x_0,y_0,z)}{q}=a^2q\left(\times 𝐀(x_0,y_0,z)\right)_z=\frac{2\pi n}{q}.$$
(6)
The crucial point is that, due to periodicity, the total flux through the system *must be zero*, i.e.

$$\underset{x,y}{}\alpha _{12}(x,y,z)=0z.$$

Consequently, the $`n`$ flux-quanta of the total flux $`2\pi n/q`$ must return in the $`+z`$ direction. This flux returns in a manner specified by the dynamics of the theory,Kajantie:1999  and it is this *response* which is the topic of interest in the current paper.

The experimental situation corresponds to applying an external magnetic field $`H`$, and then study the magnetic response of the superconductor to this field. Hence, a suitable thermodynamic description is coached in terms of a potential $`\mathrm{\Phi }(H)`$, which is a function of the *intensive* field variable $`H`$. In the simulations we have fixed $`n`$, which is analogous to fixing the magnetic induction, and a description based on the *extensive* field variable $`B`$ is more appropriate. The two approaches are related by a Legendre transformationKajantie:1999 . In principle the simulations could also be performed in an ensemble with fixed magnetic field. Technically this would be achieved by adding the term

$$HL_z2\pi n/q$$

to the action in Eq. 1. This would promote $`n`$ to a dynamical variable of the theory, and be more in accordance with the experimental situation. However, a change $`nn\pm 1`$ would require a global relaxation, and this would give very low acceptance rates, i.e. inefficient simulations.

For type-I superconductors, superconductivity vanishes for $`H>H_c`$. For type-II superconductors, a *flux line lattice* is formed at $`H=H_{c_1}`$, for smaller fields the magnetization in the sample vanishes due to the Meissner effect. By fixing $`n`$ one can not study these effect directly, however it is possible to determine a corresponding field strength from $`n`$, see Ref. Kajantie:1999, .

On the basis of simulations performed using the modified action Eq. 5, we have determined the effective temperature-renormalized interaction between two vortex-lines, and searched for the value of the GL parameter, or more precisely its square, $`x_{\mathrm{I}/\mathrm{II}}`$, where this interaction changes character from being effectively attractive to being effectively repulsive. In section II.1 we have fixed $`n=2`$ and studied the distance between the flux lines. In section II.2 we have generalized to real $`n`$, and used this to calculate a free energy difference between states containing one and two vortices. This is also a measure of the sign of the effective vortex interaction, and hence an indication of whether we have type-I or type-II behavior.

To obtain the results in Figs. 1, 2, and 5 we have performed simulations on cubic systems of size $`N=8,12,16,24,32,48`$, with $`\beta _G=1.00`$. All simulations have been performed in the broken symmetry state $`y<y_c`$, with particular emphasis on the values $`y=0.04,0.10,0.20,0.30,0.40`$. For the two largest system sizes the final datapoints are averages of approximately $`10^6`$ sweeps, whereas approximately $`10^5`$ sweeps have been performed for the four smallest system sizes.

The simulations leading to the results of Fig. 4 are quite different. They are performed for the fixed system parameters $`N=24`$ and $`\beta _G=1.00`$, and for each value of $`m`$, we have performed from $`2.510^4`$ to $`2.510^5`$ MC sweeps through the lattice. One sweep through the lattice consists of (1) conventional local Metropolis updates of $`\psi `$ and $`𝐀`$, and (2) global radial updates of $`\left|\psi \right|`$ combined with overrelaxation Kajantie:1996 ; Dimopoulos:2000  of $`\psi `$.

### II.1 Effective vortex interaction

We first clarify what is meant by effective vortex interaction in this context. In the Ginzburg-Landau model at zero temperature, one may compute a pair-potential between two vortices which consists of an attractive part due to vortex-core overlap, and a repulsive part due to circulation of supercurrents (or magnetic fields) outside the vortex core. Ignoring fluctuation effects, this furnishes an adequate way of distinguishing between type-I and type-II behavior, by asking when the attractive core-contribution dominates the magnetic field contribution or vice versa. By effective interaction, we mean a thermally renormalized pair interaction which fully takes entropic contributions into account. At low temperatures the effective interaction will revert back to the standard pair-interaction described above, but will deviate as temperature is raised, and this is particularly relevant as the critical temperature is approached, as we shall see below. We also comment further on this in the Discussion section, where we elaborate on what we perceive to be a crossover between type-I and type-II behavior.

In our simulations, the value of $`n`$ has been fixed to $`n=2`$. This corresponds to the case of two field-induced vortices which move around in the system under the influence of their mutual effective interaction. During the simulations, we have measured the transverse position $`r_{}(z)`$ of these two vortices labelled by $`1`$ and $`2`$, and the average distance between them.

$$𝖽=\frac{1}{N_z}\underset{z}{}\left|r_{}^1(z)r_{}^2(z)\right|.$$
(7)
For type-I superconductors this distance should be independent of system size, whereas for type-II we expect that this distance scales with the system size. Finite size scaling of $`𝖽`$ for various points in the $`(x,y)`$ phase diagram is shown in Fig. 1.

In the part of phase diagram which we focus on, namely the region defined by the dotted line in Fig.1 of Ref. Mo:2001, , the vortex lines are generally directed and almost straight, well defined line objects. This can be seen either by directly taking snapshot pictures of the vortex-line configurations of the system, or by computing the mean-square fluctuations around a straight-line configuration, $`|r_{}^i(z)r_{}^i(0)|^2`$, for *one* vortex line. This is in contrast to the situation in the vicinity of the critical part of $`y_c(x)`$ in Fig. 1 of Ref. Mo:2001, , where the vortex lines loose their line tension via a vortex-loop blowout Tesanovic:1999 ; Nguyen:1999 . Consequently, we can consider the theory as an effective theory for an interacting pair of straight vortex lines which interact with the dimensionless potential $`V(𝖽)`$. If we make this assumption, the probability of finding the vortices separated by a distance $`𝖽`$ in a system of $`N^3`$ lattice points with *periodic boundary conditions*, is given by

$$P_N(𝖽)=\frac{e^{V(𝖽)}\mathrm{\Omega }_N(𝖽)}{Z_𝖽},$$
(8)
where $`\mathrm{\Omega }_N(𝖽)`$ is the number of configurations with a transverse vortex-vortex distance of $`𝖽`$, and $`Z_𝖽`$ is just a normalisation factor. $`\mathrm{\Omega }_N(𝖽)`$ can be calculated, either analytically in the continuum limit

$$\mathrm{\Omega }_N(𝖽)=\{\begin{array}{cc}2\pi 𝖽\hfill & 𝖽<\frac{N}{2}\hfill \\ 2N\left(\frac{\pi }{2}2\mathrm{arccos}\left(\frac{N}{2𝖽}\right)\right)\hfill & \frac{N}{2}<𝖽<\frac{N}{\sqrt{2}},\hfill \end{array}$$
(9)
or by simple geometric counting in the case of a lattice. In the case of noninteracting vortices, i.e. $`V(𝖽)=0`$, the expectation value of $`𝖽`$ is determined only by $`\mathrm{\Omega }_N(𝖽)`$, and we find the numerical value

$$𝖽_0𝖽=\frac{1}{Z_𝖽}_0^{\frac{N}{\sqrt{2}}}𝑑𝖽\mathrm{\Omega }_N(𝖽)𝖽0.38N.$$
(10)
The separation $`𝖽_0`$ defined in Eq. 10 will be used to establish a numerical value of $`x_{\mathrm{I}/\mathrm{II}}`$. Namely, we can compute the averaged distance between vortices at fixed $`x`$ varying $`y`$, or vice versa. In the latter case, we will use the criterion that if $`𝖽`$ exceeds som value $`c𝖽_0`$ where $`c`$ is some fraction, then we have type-II behavior, otherwise it is type-I. The quantity $`𝖽`$ at fixed $`y`$ will turn out to be an S-shaped curve as a function of $`x`$, increasing from small values to large values as $`x`$ is increased. We interpret this as yet another manifestation of the crossover from type-I to type-II behavior, and we have chosen to locate the crossover region $`x_{\mathrm{I}/\mathrm{II}}`$ at the value of $`x`$ where the curves change most rapidly, which is roughly when $`𝖽𝖽_0/2`$. As we shall see (see Fig. 6), different crossover criteria give consistent results. The quantity $`P_N(𝖽)`$ can be estimated from histograms, and then we can use Eq. 8 to determine the pair potential. Depending on whether we consider type-I or type-II superconductors we expect to see an attractive or a repulsive potential. Fig. 2 shows the potential $`V(𝖽)`$ for the same points of the phase diagram as Fig. 1.

### II.2 Free energy

In Eq. 5 we have used $`n`$ to indicate an integer number of flux tubes, but in principle there is no reason to limit $`n`$ to integer values, and we will use $`S[\alpha ,\psi ,m]`$ to denote a generalisation to real $`n`$.

We have considered the free energy difference between a state containing zero vortices, i.e. $`m=0`$ and a state containing $`n`$ vortices. We can not measure absolute values of the free energy, but by differentiatingKajantie:1999

$$e^{F(m)}=\mathrm{𝖳𝗋}e^{S(m)}$$
(11)
with respect to $`m`$, and then integrating up to $`n`$, we can calculate $`\mathrm{\Delta }F(n)=F(n)F(0)`$,

$$\begin{array}{cc}& \frac{\mathrm{\Delta }F(n)}{L_zq^2}=\hfill \\ & 2\pi \beta _G_0^n𝑑m\left[\underset{W(m)}{\underset{}{2m+\frac{1}{\pi N_z}\underset{z}{}\alpha _{12}(x_0,y_0,z)_m}}\right].\hfill \end{array}$$
(12)
To calculate $`\mathrm{\Delta }F`$, we have then varied $`m`$ in steps of $`\mathrm{\Delta }m=0.05`$, and performed the integration in Eq. 12 numerically. Using shift symmetries, it can be shown Kajantie:1999  that $`W(n)`$ is equal to

$$W(n)=\frac{2n}{N_xN_y},$$
(13)
and the behavior for intermediate real values is shown in Fig. 3. Increasing $`n`$ from $`0`$ to $`1`$ costs a free energy $`\mathrm{\Delta }F(1)`$, and adding two vortices costs an amount $`\mathrm{\Delta }F(2)`$. We will *always* have $`\mathrm{\Delta }F(2)>\mathrm{\Delta }F(1)`$, but the question is whether $`\mathrm{\Delta }F(2)2\mathrm{\Delta }F(1)`$. We may regard $`F(n+2)+F(n)2F(n+1)`$ as the discrete second derivative of the free energy with respect to particle number, which is nothing but the inverse compressibility $`K^1`$ of the vortex-system. In the thermodynamic limit this quantity can never become negative. However, its vanishing signals the onset of *phase-separation* of the vortex system, which we again interpret as a lack of stability of the vortex-loop tangle, characteristic of type-I behavior.

### II.3 Vortex compressibility, separation, and crossover

We next define a quantity $`\mathrm{\Delta }T`$ by the relation

$$\mathrm{\Delta }T=\frac{\mathrm{\Delta }F(1)}{L_z}\frac{\mathrm{\Delta }F(2)}{2L_z},$$
(14)
which means that $`\mathrm{\Delta }T`$ measures the relative free-energy difference between adding one vortex to the system and half of that adding two vortices to the system. Intuitively it is therefore clear that it measures the sign of the vortex interactions, and hence determines whether we are in the type-I or type-II regime. $`\mathrm{\Delta }T>0`$ signals attractive interactions, i.e. type-I behavior, whereas $`\mathrm{\Delta }T0`$ signals repulsive interactions, i.e. type-II behavior. We have calculated $`\mathrm{\Delta }T`$ by using Eq. 12 and Eq. 14, the results are shown in Fig. 4. The main qualitative result from these simulations is again that $`x_{\mathrm{I}/\mathrm{II}}(y)`$ is a declining function of $`y`$. Note also that $`K^1=\mathrm{\Delta }T`$, and hence a positive $`\mathrm{\Delta }T`$ clearly implies phase-separation and instability of the vortex system, characteristic of type-I behavior. This is precisely what we see for small $`x`$ in Fig. 4.

Finite size scaling of $`𝖽`$ and studies of $`\mathrm{\Delta }T`$ differentiate nicely between strongly type-I and type-II superconductors, but it is difficult to locate a value of $`x_{\mathrm{I}/\mathrm{II}}(y)`$ with any great precision. Fig. 5 shows $`𝖽(x)`$ for different values of $`y`$, along with a horizontal line at $`𝖽_0/2`$, where $`𝖽_0`$ is the average separation between vortices had they been non-interacting. We have found that $`𝖽_00.38N`$ in our simulations. We have, rather arbitrarily, taken the interception of this horizontal line with the curve $`𝖽(x)`$ as $`x_{\mathrm{I}/\mathrm{II}}`$.

The curves of $`𝖽(x)`$ do not get significantly sharper with increasing system size, and there are no particular sharp features in $`S[\alpha ,\psi ,2]`$ as $`x`$ is increased beyond $`x_{\mathrm{I}/\mathrm{II}}`$. Fig. 6 shows the intercepts from Fig. 5. Due to the features in the curves of Fig. 5, and how the results of Fig. 6 are obtained from them, we tentatively conclude that the computed line of Fig. 6, corresponding to the dashed line of Fig. 1 in Ref. Mo:2001, , is a *crossover* and *not* a phase transition. However, we comment further on this in the concluding section.

As already indicated, there is some arbitrariness in the location of $`x_{\mathrm{I}/\mathrm{II}}(y)`$ in Fig. 6, however the four points labelled by (IA, IB) and (IIA, IIB) clearly are in the type-I and type-II regimes, respectivly. This is demonstrated in Figs. 1 and 2.

## III Discussion

From Figs. 1 and 2 we conclude that there is a crossover line separating effective attractive vortex interactions from effective repulsive ones, i.e. type-I and type-II. This line can either be crossed by changing $`x`$, i.e. IA $``$ IIA, or by changing the temperature i.e. IB $``$ IIB in Fig. 6. This means that for $`x`$ values in a suitable range, we can have in principle have a *temperature induced* crossover from type-I to type-II superconductivity. Finally, we note that $`x_{\mathrm{I}/\mathrm{II}}(y)`$ deviates significantly from the mean-field value of $`x_{\mathrm{I}/\mathrm{II}}=0.5`$.

Deep in the type-I regime, we find clear evidence of *attractive* interactions. In the type-II regime the repulsive interactions appear to be weak, and the results are essentially also consistent with two *randomly placed* vortex lines, i.e. not interacting, but not consistent with an attractive force between the vortices. Therefore, what the results unequivocally show is that by fixing the material parameter $`x`$ and varying the temperature-like variable $`y`$, the character of the effective pair-potential is altered inside the superconducting regime, significantly away from the critical line.

### III.1 Simulations

The simulations with $`n=2`$, and the simulations with a real $`m[0,2]`$ are quite different, and we will discuss them in turn. An important feature of *all the simulations* in the present work is *slow dynamics*.

For the $`n=2`$ simulation, where we have monitored $`𝖽`$, we find that deep in the type-I regime the simulations are quite straightforward, the vortices stay close together with only small fluctuations, and a moderate number of sweeps is sufficient to get good statistics. However when we increase $`x`$ towards $`x_{\mathrm{I}/\mathrm{II}}`$ the effect is *not* that $`𝖽`$ stabilises at a higher value, instead we get get fluctuations between a type-I like state where the two vortices are close together, and a type-II like state with large vortex-vortex separation. This picture persists as $`x`$ is increased into the type-II regime, the only difference is that the fraction of time spent in the type-I like state state decreases. In fact the results of Fig. 1 and Fig. 2 are in the type-II region quite close to what we would get from noninteracting vortices. What makes the simulations difficult is that moving a vortex line in the transverse direction is a *global* change, and thereby very slow. Timeseries of $`𝖽`$ show characteristic time scales of $`10^4`$ sweeps, so long simulations $`10^6`$ sweeps over the lattice are required to obtain acceptable accuracy. A truly high-precision determination of $`x_{\mathrm{I}/\mathrm{II}}(y)`$ would surely benefit from a specialized algorithm for the MC updates.

To get good results one should take the $`N\mathrm{}`$ limit. The conclusions from the results of Fig. 1 are based on this limit, whereas those drawn from Figs. 2,5 and 6 are based on the fixed system size $`N=48`$. We have not performed a systematic study of finite size effects, but the curves in Fig. 5 *do have* subject to finite size effects in them. The trend is that curves move to the right upon increasing system size, this very likely explains the apparent discrepancy between the tricritical point (where the $`N\mathrm{}`$ limit has been applied), and the remaining points in Fig. 6.

Note that Eq. 1 is a continuum field theory, and $`x_{\mathrm{I}/\mathrm{II}}(y)`$ is *not* a critical point, hence the continuum limit $`\beta _G\mathrm{}`$ should be taken. Our experience from the large-scale simulations performed in Ref. Mo:2001,  indicates that $`\beta _G=1.00`$ provides conditions in the simulations already quite close to the continuum limit. We have therefore chosen to work with $`\beta _G=1.00`$ and focused our efforts on considering large systems and long simulations.

The Monte-Carlo computations of $`\mathrm{\Delta }T`$ have been even more time consuming, because we have had to do the simulations for 41 different values of $`m`$. We have therefore limited ourselves to considering only the system $`L=24,\beta _G=1.00`$, for a discussion of finite $`N`$ / finite $`\beta _G`$ effects see Ref. Kajantie:1999, . The relaxation time for these simulations has been particularly long in the limits $`m1^{}`$ and $`m2^1`$, and we therefore have performed much longer simulations in these limits than for intermediate values of $`m`$. From Fig. 3 and Eq. 13 it is seen that $`W(m)`$ is a quantity of order $`𝒪(N^2)`$ whose finiteness originates in the difference between two $`𝒪(1)`$ quantities. Consequently, it is difficult to get numerically precise results. This has in particular been the case in the limits $`m1^{}`$ and $`m2^{}`$.

### III.2 Crossover

The physical picture that emerges for the thermal renormalization of the vortex interactions, is the following. At low temperatures, for the $`\kappa `$-values we consider, the system is in the type-I region with a fairly deep minimum in the effective pair-potential between vortices at short distances, leading to an attractive interaction. At very short distances, we find a steric repulsion on the scale of the lattice constants in the system due to the large Coulomb barrier that must be overcome to occupy a link with two or more elementary vortex segments. This lenght scale represents the size of the vortex core in the problem. Upon increasing the temperature to the vicinity of the line $`y_c(x)`$, we do not find large transverse meanderings of the individual vortex lines as we move along each vortex line, rather the vortex lines are essentially straight. Therefore, we believe that it is not entropic repulsion due to the bare steric repulsion in the problem, of the type which it seems reasonable to invoke for strongly fluctuating elastic strings Zaanen:2000 ; Mukhin:2001  that renormalizes the vortex interactions in the way that is seen in Fig. 2. Rather, what appears to happen is that the vortex lines slosh back and forth in the minimum of the effective potential well as essentially straight lines. Hence, to a larger and larger extent as temperature is increased, they experience the hard wall in the interactions at small distances, and the weak attraction at large distances, effectively washing out the minimum in the potential, thus making it effectively more repulsive.

This is also seen in our simulations (not shown in any of the figures) when we monitor the transverse meandering fluctuations of each vortex line, $`|r_{}(z)r_{}(0)|^2`$, as well as the mean square fluctuations of the intervortex distance, $`𝖽^2𝖽^2`$, where $`𝖽`$ is defined in Eq. (7). The former is small deep in the superconducting regime, and remains small as the line $`y_c(x)`$ in Fig. 1 of Ref. Mo:2001,  is approached, while the latter increases dramatically as the dotted crossover line is crossed. It is precisely this fact which makes the simulations extremely time consuming.

One should however keep in mind that, since we are considering the full GL theory in our simulations, and not the linearized London limit, it is in any case a drastic simplification to view the effective interaction between vortices as a simple pair potential.

Finally, we note that, although our present simulations, which by necessity are on finite-sized systems, indicate that the change from type-I to type-II behavior is a crossover, we cannot rule out the possibility that it is elevated to a true phase-transition in the thermodynamic limit. More work is needed to clarify if this is indeed the case, but this will have to await the next generation of massive parallell computers. Questions that need to be addressed in this context, are: What is the order parameter of such a transition, and what symmetry, if any, is being broken.

## IV Conclusion

We have considered the effective interaction between two vortices in the full GL model, and how this effective interaction is influenced by thermal fluctuations. We have included fluctuations in the gauge fields, as well as the phase- and amplitude-fluctuations of the complex scalar matter field of the problem. We have found that the effective interaction changes from being attractive to being repulsive at $`x_{\mathrm{I}/\mathrm{II}}`$. This means a change from type-I to type-II behavior. We have found that $`x_{\mathrm{I}/\mathrm{II}}`$ is below the standard quoted value of $`0.5`$, and is a function of the temperature-like parameter $`y`$. This means that at the critical point, the value of the GL parameter that separates type-I from type-I behavior is smaller than $`1/\sqrt{2}`$. The line $`x_{\mathrm{I}/\mathrm{II}}(y)`$ appears to be a crossover, and not a true phase transition. The above seems to offer a simple explanation for the experimental observation that elemental Ta and Nb superconductors show a crossover from type-I to type-II behavior as the temperature is increased towards $`T_c`$. Previous explanations based on mean-field theories and not involving thermal fluctuations required two additional temperature dependent $`\kappa `$-values to be defined Maki:1969 .

We acknowledge support from the Norwegian Research Council via the High Performance Computing Program (S.M.,J.H.,A.S.), and Grant Nos. 124106/410 (S.M.,A.S.) and 148825/432 (A.S.), and A.S. thanks E. H. Brandt for useful comments. J.H. also acknowledges support from NTNU via a university fellowship.