CHAPTER 1 INTRODUCTION
1.4 Goals of the present work
Although most people believe the existence of the Abrikosov lattice phase that solved from GL equation by LLL approximation, but T. Giamarchi and P. Le Doussal disagree([19] T. Giamarchi). They claimed there’s no Abrikosov phase by two approaching: gauge glass model([20]M. P. A. Fisher; [21] D. S. Fisher) and elastic lattice structure at small scale([22]M. Feigelman). Fig. 1-4-1 is the phase diagram of their prediction. I tried to plot the phase diagram by using Monte Carlo simulation. I analyzed the results of the MC simulation directly in mathematical way to check the phase diagram.
- 8 -
Fig. 1-4-1: The phase diagram by T. Giamarchi and P. Le Doussal’s prediction([19] T. Giamarchi).
Comparing with Fig. 1-3-1, one can find there’s no Abrikosov phase in this phase diagram.
- 9 -
Chapter 2
Description of a layered superconductors in
strong magnetic fields by the LLL approximation
2.1 Landau levels and the quasimomentum basis
We focuses on the properties of the two dimensional type II superconductors with the external magnetic field in the region H𝑐1≪ H < H𝑐2 in the thesis. LLL approximation is the most acceptable theory where H ≈ H𝑐2. In LLL theory, since magnetization is small we replace the field inside superconductor B by external field H which is essentially homogeneous. So we can drop the nonlinear term of the Ginzburg-Landau equation.
As the thermal fluctuations becoming strong enough, excitations of the lattice are no longer invariant under the symmetry transformations. So we need to regauge the system with LLL theory by the quasimomentum basis.
Solve the linear Ginzburg-Landau equation order parameter expanding in quasimomentum basis within the LLL([23] D. Li and B. Rosenstein)
ψ(x, y) = ∑ 𝐶𝑘𝜑𝑘
- 10 -
where the coefficients 𝐶𝑘 are complex numbers and 𝜑𝑘 are quasi-momentum basis. We can find the property of 𝜑𝑘
𝜑𝑘 = exp *−𝑖𝑥𝑘𝑥+𝜑0(𝑥 − 𝑘𝑦, 𝑦 + 𝑘𝑥). (2.1.3) It is proved that quasi-momentum basis satisfy magnetic transitions([24] B.
Rosenstein) which is defined as
T𝑑𝜑𝑘= 𝑒𝑖𝒌𝒅𝜑𝑘, (2.1.4)
here 𝐝 is the general displacement vector and T𝑑 is magnetic translation operator.
The sample size in the simulation is finite and had following dimensions([25] H. Y.
Lin):
We work in reciprocal lattice vector, so that 𝒌 with the basis vector is
𝒌 = k1𝒅̃𝟏+ k2𝒅̃𝟐 (2.1.8)
Here k1 = 0,1, … ,L 1 and k2 = 0,1, … ,L 1 due to we choose the lowest Landau level wave function with quasi-momentum 𝒌.
- 11 -
Fig. 2-1-1: Interpretation of the lattice vector 𝑑1 and 𝑑2. 𝑎 is the distance of a vortex to another adjacent vortex.
The quasi-momentum basis satisfy magnetic translations
T 𝜓(𝑥, 𝑦) = 𝑒𝑖𝑘 𝜓(𝑥, 𝑦), (2.1.9)
and then we have
𝑒𝑖𝑘 = 1
k𝑥 =2𝜋
𝐿 𝑛𝑥, 𝑛𝑥 = 0, ±1, ±2, …
k𝑦 =2𝜋
𝐿 𝑛𝑦, 𝑛𝑦 = 0, ±1, ±2, … (2.1.10)
Because of using the LLL approximation, the ranges of 𝑛𝑥 and 𝑛𝑦 are 𝑛𝑥 = 0,1, … , 𝐿 − 1
𝑛𝑦 = 0,1, … , 𝐿 − 1 (2.1.11)
- 12 -
2.2 Free energy
We start at the two dimensional GL free energy:
F = ∫ 𝑑𝑥𝑑𝑦 ℏ where 𝜓 is the order parameter of the superconductivity, 𝑨 is the magmetic vector potential. 𝑚 and e∗ are the mass and charge of the Cooper pair. 𝛼′(𝑇) = 𝛼𝑇𝑐(1 − 𝑡) and 𝑏′(𝑇) are phenomenological parameters. 𝑈(𝑥, 𝑦), the disorder term, is like ( ) which is the disorder coefficient by Larkin ([18] A. I. Larkin).
𝑈(𝑥, 𝑦)𝑈(0,0)
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ = 𝑅𝛿(𝑥)𝛿(𝑦) (2.2.2)
Let us focus the clean system first. The GL free energy equation for pure vortex system is
In order to get the dimensionless LLL free energy, we need to rescale ψ2 →
√𝑙2𝑏′4𝜋𝑇 ψ2, x → x𝑙, y →y𝑙 and 𝑙 = ξ
Here is the GL free energy in the form that we'll use it in Monte Carlo simulation.
The Gaussian integrals referred to the paper of B. Rosenstein and D. Li ([24] B.
- 13 - with decomposition of arbitrary momentum 𝑞 into its rational part k, which belongs to the Brillion zone, and an integer part 𝑸, belonging to the reciprocal lattice
𝒒 = 𝒌 + 𝑸, 𝑸 = 𝑄1𝒅̃𝟏+ 𝑄2𝒅̃𝟐 (2.2.8)
The inverse Fourier transform of Eq. (2.2.6) is 𝜑(𝒓)𝜑𝒌∗(𝒓) = ∑ exp,𝑖(𝒌 + 𝑸) ⋅ 𝒓- exp [𝜋𝑖
Since the superfluid density and its Fourier transform defined as:
𝜌(x, y) = |𝜓(𝑥, 𝑦)|2 = ∑ 𝐶𝑘𝐶𝑙𝜑𝑘𝜑𝑙
𝑘,𝑙
(2.2.10)
𝜌̃(𝐩 + 𝐏) = 1
2𝜋𝐿2∫ 𝑒𝑥𝑝,−𝑖(𝒑 + 𝑷) ⋅ 𝒓- × 𝜌(𝑥, 𝑦) 𝑑𝑥𝑑𝑦 (2.2.11) From the equations from Eq. (2.2.6) to Eq. (2.2.11) we get the quadratic term in Eq.
(2.2.4)
while the quartic term takes a form
- 14 -
The disordered term of the GL free energy is
∫ 𝛼𝑇𝑐(1 − 𝑡)𝑈(𝑥, 𝑦)|𝜓|2
The detail definition of Eq. (2.2.15) to Eq. (2.2.15) are shown in Appendix A.
- 15 -
2.3 Disorder parameter 𝜻 of Li and Nattermann
M. S. Li and T. Nattermann([26] T. Natterman) defined the dimensionless parameter ζ to control the relative disorder strength, and expand the random Gaussian disorder in renormalized Hermite polynomials to express the disorder term of GL free energy equation. The disorder term of M. S. Li and T. Nattermann is
∫ 𝑑2𝑟 𝛼𝛿𝑇𝑐(𝑟)|𝜓|2 (2.3.1) They defined the disorder parameter with following relation
ζ̃ = 𝜁𝑏12 𝜋12(1 − 𝑡 − 𝑏)
(2.3.5)
Connecting our disorder term with their notation, we have 𝑅𝑙𝛼𝑇𝑐(1 − 𝑡)
8𝜋𝐿2(𝜋𝑏′𝑇)12
= 𝑙2𝑎𝑇2𝜁̃2
64𝜋2𝐿2 (2.3.6)
By controlling the specific ζ̃ to generate the corresponding complex random fields, we can get various degree disorder vortex systems.
- 16 -
2.4 MC simulation method.
Metropolis algorithm:
The standard Monte Carlo method with Metropolis algorithm ([28] D. P. Landau and K. Binder; [27] Y. Kato) was used to simulate the two-dimensional pure and disordered vortex system. In the classic Metropolis method, we use a transition probability which depends on the difference of energy ΔE between the initial and trial configuration to determine whether the trial configuration is accepted or not.
Now I introduce the Monte Carlo method as follows. First, we choose an initial configuration and calculate the initial energy E𝑚. Second, we choose a site C𝑗 ∈ 𝐶𝑁 randomly and generate the trial configuration with C𝑗𝑛𝑒𝑤 by using the rule:
C𝑗𝑛𝑒𝑤 → C𝑗𝑜𝑙𝑑+ ϵΔC, where ΔC is a complex number which is chosen randomly from the region |ReΔC| ≤ 1 and |ImΔC| ≤ 1 in the complex plane. Third, we calculate the energy E𝑛 of trial configuration and the difference of energy ΔE, here
ΔE = E𝑛− E𝑚. If ΔE ≤ 0, the system accepted the trial configuration, but if ΔE ≥ 0, the trial configuration is accepted with a probability exp(−βE). Generating a
random number uniformly in the interval ,0,1-, if ≤ exp(−βE) the trial configuration is accepted, otherwise it is rejected. This process is called Monte Carlo step/site (MCS/site). Note that the old configuration is still counted again for
averaging if the trial configuration is rejected. By using Monte Carlo method, the system will fall into the stable states and reach the equilibration, and the
characteristics of vortex system can be measured.
- 17 -
Monte Carlo calculations:
We used
C𝑙𝑛𝑒𝑤1,𝑙2 = C𝑙𝑜𝑙𝑑1,𝑙2+ 𝛿𝑙1 𝑗1𝛿𝑙2 𝑗2Δ (2.4.1) to vary the value of a specific wave function coefficient C𝑗 in our Monte Carlo simulation, here Δ = ϵΔC. Note that j and l are vectors which composed of two reciprocal vectors d̃₁ and d̃₂. Furthermore we used these equation to calculate the energy of trial configuration and only discussed the changes of the summation of wave function coefficient product, the detail of Monte Carlo calculations are worked out in Appendix C. The summation of wave function coefficient product of trial configuration is hence we can store it to simulate the vortex system more efficiently. The old
calculation results always can be applied in new one and a lot of computer time is saved, the CPU time in one Monte Carlo step ∝ L² . There are 16 × 16 numbers of vortices in our simulation. We took 6 × 10⁶ MC steps to reach the thermal
equilibration and calculated the physical quantities over 1 × 106~1 × 10⁷ MC steps.
- 18 -
The physical quantities were measured every 30 ~50 MC steps. We control ϵ in a reasonable region to make the acceptance ratio is 0.3 ~0.4 and then the vortex system reach the thermal equilibrium state efficiently. All the simulations were started from the heating processes with the initial configuration which is defined as
follows: C𝑖 = √| 𝛽𝑇|
𝐴, here C𝑖 is one of all coefficients of wave function and others are equal to zero, 𝛽̅𝐴 ~ 1.16 is the mean-field value of the Abrikosove ratio.
- 19 -
Chapter 3
Indicators of the four vortex phases and the MC data Analysis
3.1 Vortices configuration
Fig. 3-1-1 are the snapshots of the spatial distribution of the order-parameter field |𝜓(𝑥, 𝑦)|2 for 𝑎𝑇 < 𝑡𝑚 and 𝑎𝑇 > 𝑡𝑚, respectively. There are 16 × 16 = 256 vortices in each sample and we used 6 × 106 Monte Carlo steps to make sure the system reaches the equilibration by referring the MC simulation by Y. Kato and N.
Nagaosa ([27] Y. Kato). The melting temperature was calculated of the sample with finite size,
𝑡𝑚 = −13.02 for N𝑠 = 256 (3.1.1)
Fig. 3-1-1(a): snapshots of the spatial distribution of the order-parameter field
|𝜓(𝑥, 𝑦)|2 for 𝑎𝑇 < 𝑡𝑚.
- 20 -
Fig. 3-1-1(b): snapshots of the spatial distribution of the order-parameter field |𝜓(𝑥, 𝑦)|2 for 𝑎𝑇 > 𝑡𝑚.
One can easily see that the vortices arrayed regular and randomly for 𝑎𝑇 < 𝑡𝑚 and 𝑎𝑇 > 𝑡𝑚 in Fig. 3-1-1. The Fourier transform of Fig. 3-1-1 is shown in Fig. 3-1-2. One can find a big change while the pattern is structure less above 𝑡𝑚 and the pattern below 𝑡𝑚 shows the six sharp peaks with hexagonal symmetry indicating the existence of the lattice of the vortices. I try to analysis the change of the Fourier transform pattern for disorder parameter ζ and reduced temperature 𝑎𝑇 due to ζ = 0~1 and 𝑎𝑇 = −17~ − 11 to find the phase transition.
- 21 -
Fig. 3-1-2(a): The Fourier transform of Fig. 3-1-1(a).
One can see beautiful hexagonal lattice
Fig. 3-1-2(b): The Fourier transform of Fig. 3-1-1(b).
It’s a structure less pattern instead of beautiful hexagonal lattice
- 22 -
3.2 Rotation average
We made discrete data with reduced temperature 𝑎𝑇 = −17~ − 11 and reduced disorder ζ = 0~1 which are snapshots in the equilibrium state. There are 20 samples with each 𝑎𝑇 and ζ. The first thing we need to do is averaging the 20 samples of the Fourier transform pattern. Since they are just snapshots in
equilibrium state, averaging them let us get rid of the extreme un-objective sample so that we can mark the phases correctly. Consider the patterns with beautiful lattice for 𝑎𝑇 < 𝑡𝑚 first. I rotated the highest peak of the first ring to the same angle and averaged them so the information of the patterns will not lose in averaging directly.
It’s no doubt that the other peaks will be moved to the same coordinate because it’s always a hexagonal lattice at first ring. And I averaged the six peaks to eliminate the height difference between the peaks because I don’t want to always see a peak much taller than others due to I made it in the rotation. Fig. 3-2-1 is the simple
interpretation of the rotation. Comparing the patterns before and after rotation averaging, one should find they are the same since it is structure less originally and they have a big difference since it is a hexagonal lattice originally. Fig. 3-2-2 shows the series procedure of the rotation averaging.
- 23 -
Fig. 3-2-1: Interpretation of the rotation: I rotated all the peaks to the same position.
Fig. 3-2-2(a): The pattern of vortices averaged directly. One will get averaged patterns like this with any reduced temperature and disorder without rotation.
- 24 -
Fig. 3-2-2(b): There is a tallest peak much taller than others due to that I choose the highest peak rotating to the same position.
Fig. 3-2-2(c): After angle averaging: averaged six parts of dividing by every 60 degrees to eliminate the un-objective peak.
- 25 -
I concluded the trend of the patterns changing. As the temperature increasing, the vorteice melted by thermal fluctuation. And as the disorder increasing, the vortices approach pinning on the disorder and also destruct the lattice. Fig. 3-2-3 and Fig. 3-2-4 show the changing with reduced temperature and reduced disorder.
Fig. 3-2-3: Series Fourier transform pattern, the temperature getting higher from (a) to (d).
- 26 -
Fig. 3-2-4: Series Fourier transform pattern, the disorder getting higher from (a) to (d).
According to the appearance of the samples after averaging, I marked the phase diagram to three parts: the Fourier transform pattern with (a)beautiful hexagonal lattice, (b)structure less lattice and (c)disturbed but still looked like hexagonal lattice that show on Fig. 3-2-5. But it’s difficult to separate the phases precisely so I tried some numerical analysis that I introduce in the following section.
- 27 -
Fig. 3-2-5: The phase diagram can be separate to three parts roughly. In solid phase, there are six sharp peaks separately formed the beautiful hexagonal lattice. While in vortex liquid phase, there’s nearly no peak observed in the diagram so looked like a smooth ring. And one can see the vortices “sticky” connected by each other in the glass phase.
3.3 Indicator of melting line
To determine the melting line precisely, we need a numerical analysis. I integrated the first ring and the six peaks individually as shown in Fig. 3-3-1. By comparing the integration of each individual peaks and the first ring
∫
𝑖𝑛𝑡𝑒𝑔𝑟𝑎𝑙,𝑝𝑒𝑎𝑘𝑠-∫ 𝑖𝑛𝑡𝑒𝑔𝑟𝑎𝑙,𝑟𝑖𝑛𝑔- ≡ 𝑅𝑚 (3.3.1)
- 28 -
Fig. 3-3-1: Interpretation of R: (a) is enlarged pictures of a single peaks, (b) is the snapshot including the first ring only.
In theoretical prediction, there's almost no difference between these two integral in solid. So I expected there's a big step at the melting line. First I choose the highest point in the first ring, and summed over the value of nearby points. The other five peaks are determined by choosing the five peaks having a included angle respond to origin and the highest point in 60, 120, 180, 260 and 320 degrees. For ideal solid phase the ratio 𝑅𝑚 must be 1. And the ratio for ideal liquid is about 0.42 that I tested it by a simulation of ideal liquid as Fig. 3-3-2.
- 29 -
Fig. 3-3-2: Ideal liquid Fourier transform of |𝜓|2
3.4 Indicator of glass (irreversibility) line
To determine the glass line precisely, we need to compare the vortices position in two different initial conditions. If the system is in Brrag glass phase or vortex glass phase. The vortices will prefer pinning on the disorder. So I rotated the highest votex peak of the first ring to the same position with two different initial conditions like I did in the average. And I multiplied these two plots and integrated it. Before
integrating, I had eliminated the points value nearby the highest point included the highest point to avoid always getting high integral value nearby the highest point. If the system isn't in glass phases the peaks with two initial conditions will be in almost the same position, and the integral must be much larger than the integral with
- 30 -
system in glass phases. The precede work before the integration show in Fig. 3-4-1.
For convenient, I call the integration as 𝐼𝑔.
Fig. 3-4-1: The precede work before the integration: (a) vortices with random distribution disorder; (b) vortices with lattice distribution disorder;
(c) after multiplied of (a) and (b); (d) the highest peak is eliminated.
- 31 -
Chapter 4
The phase diagram of the vortex matter
4.1 Melting line
Using the indicator defined in section 3.3, I sketched the 𝑎𝑇− R𝑚 diagram of the pure system in two initial conditions (In ideal cases, there’s no difference between two systems with random distribution of disorder and lattice like
distribution of disorder.)(Fig. 4-1-1). From Fig. 4-1-1 we can see that the boundary value of solid and liquid phase is about 0.5. If the value of R𝑚 is above 0.5, I said the system is more like in solid phase. Otherwise, if R𝑚 is below 0.5, the system is more like in liquid phase.
Fig. 4-1-1(a): 𝑎𝑇− R𝑚 diagram with random distribution of disorder in pure system.
m
- 32 -
Fig. 4-1-1(b): 𝑎𝑇− R𝑚 diagram with lattice distribution of disorder in pure system.
Use R𝑚 = 0.5 to be the indicator of the boudary of solid phase and liquid phase. I calculated the ratio R𝑚 of each samples with reduced temperature during 𝑎𝑇 = −17~ − 11 and disorder parameter during ζ = 0~1. As shown in Fig. 4-1-2, it's clear to see the melting line in the 𝑎𝑇− ζ diagram. It matches the melting line of the contour plot of 𝑎𝑇− ζ diagram that one can find the transition line at the same position (Fig. 4-1-3).
Fig. 4-1-2(a): If the ratio R of the system is larger than 0.5, it marked red point. In another way, if the ratio R of the system is smaller than 0.5, it marked blue point. It's the statistical of R with random distribution initial condition.
m
Rm
aT
- 33 -
Fig. 4-1-2(b): It's the statistical of R with lattice-like distribution initial condition. One can see the same melting line clearly both in random distribution i.c. or in lattice-like distribution i.c..
Fig. 4-1-3(a): Contour plot of 𝑎𝑇− ζ diagram of R𝑚 with random distribution disorder
aT Rm
aT Rm
- 34 -
Fig. 4-1-3(b): Contour plot of 𝑎𝑇− ζ diagram of R𝑚 with lattice like distribution disorder. one can see the transition line is the same with two different distribution of disorder by comparing (a) and (b), and fitted the transition line in Fig. 4-1-2.
4.2 Glass (irreversibility) line
Using the indicator defined in section 3.4 and sketched the 𝑎𝑇− 𝐼𝑔 diagram directly without considering the thermal fluctuation is strong at high temperature, we got a glass line isn't objective (Fig. 4-2-1) that the glass point should be higher as the temperature getting higher because the peak is hard to pin on the disorder by the strong thermal fluctuation at high temperature. The reason to see a un-objective trend is because the vortices almost melted at high temperature, so of course the integral 𝐼𝑔 is smaller than the vortices integral under low temperature (no sharp
Rm
aT
- 35 -
peaks).
Fig. 4-2-1: Density plot without correction, one can see the integration I getting higher as temperature getting higher at low temperature but getting lower as temperature getting up.
Here are two ways to solve the problem. One is see the plot fixed one
temperature at a time. The border on the step is the transition point of the glass line (Fig. 4-2-2). By observing the plot fixed one temperature at a time, one will not confused by the number but focuses on the variation of 𝐼𝑔. Another way to see the trend of glass line is dividing the integral value 𝐼𝑔 of disordered system by the integral 𝐼𝑔𝑝𝑢𝑟𝑒 of pure system. In ideal case, vortices with two different initial
condition in pure system map to each other perfectly and the integral is always larger
aT Ig
- 36 -
than the case in disorder system. So it's reasonable to divided the integral value of disorder by the integral of pure system to "normalize" the integral values. We can see the trend in the contour plot after correction by the second way in Fig. 4-2-3 that it fit the transition points in Fig. 4-2-3 with different 𝑎𝑇.
Fig. 4-2-2: ζ − 𝐼𝑔 diagram with fixed 𝑎𝑇 at a time. The decrease speed as disorder getting larger of 𝐼𝑔 is slow down at high temperature. It means it's hard to pin on the disorder at high temperature so there's no big difference between low ζ and high ζ.
Ig
ζ
- 37 -
Fig. 4-2-3: Contour plot after correction. One can see the integration 𝐼𝑔 getting higher sharply at 𝑎𝑇 = −13.
4.3 Phase diagram
Assemble the results in section 4.1 and 4.2. I got a 𝑎𝑇− 𝐼𝑔 phase diagram roughly showed in Fig. 4-3-1. The melting point of my analysis in pure system is at 𝑎𝑇 = −13 approximately. It’s the same as the calculation by Kato ([27] Y. Kato) in MC simulation with finite sample. The trend of melting line as disorder getting large fit the theory prediction that the melting point with larger disorder is at lower
temperature. Because the well accuracy of melting line, the glass line in my analysis is reliable. From my analysis, the Abrikosov lattice is existence. There is no pinning effect since the disorder is weak. From Fig. 4-2-2 one can see there is a big step when
ζ
aT
- 38 -
the disorder get strong enough. Also from the big difference of patterns of the Fourier transform of |𝜓|2 between solid phase (Abrikosov phase) and glass phase (Brag glass or vortex glass) we saw in Fig. 3-2-4, the Abrikosov lattice is existence indeed.
Fig. 4-3-1: 𝑎𝑇− 𝜁 phase diagram. The blue line is the glass transition line; The red line is the melting transition line.
- 39 -
Chapter 5
Conclusion
The vortex simulation in highly anisotropic layered type II superconductors has been studied by Monte Carlo simulation in two dimensional Ginzburg-Landau model with the quasi-momentum basis. Vortices structure are studied with disorder
parameter ζ = 0~1 and reduced temperature 𝑎𝑇 = −17~ − 11 in the thesis. I developed the rotation averaging to analysis the snapshots of the Fourier transform of the superfluid density. Using the rotation averaging, I compared the diagrams of the average of samples with each disorder parameter and reduced temperature before and after rotation to classify the diagrams into three categories: Abrikosov lattice phases, vortex liquid phase and glass phase. In solid phase (Abrikosov phase), there are six sharp peaks separately formed the beautiful hexagonal lattice. While in vortex liquid phase, there’s nearly no peak observed in the diagram so looked like a smooth ring. And one can see the vortices “sticky” connected by each other in the glass phase. To identify the transition line precisely, I made two indicators for melting line and glass line as R𝑚 and 𝐼𝑔 to analysis phase diagram numerically. R𝑚 is defined by the compare of the integral of finite area near ideal solid peaks and the integral of the first ring (the central lattice). R𝑚 must be 1 in ideal solid phase since all contribution of the integral of the first ring comes from the six sharp peaks only. In another way, R𝑚 must be small because most contribution of the integral of the first ring doesn’t come from the area near the peaks even there’s no peaks in the
- 40 -
liquid phase. The glass transition line indicator 𝐼𝑔 is defined by the multiplication of the samples at the same disorder parameter and reduced temperature but with two
liquid phase. The glass transition line indicator 𝐼𝑔 is defined by the multiplication of the samples at the same disorder parameter and reduced temperature but with two