# Vortex Lattice Melting Line in Superconductors with Paramagnetic Pair-
Breaking

Dai Nakashima and Ryusuke Ikeda  Department of Physics, Kyoto University,
Kyoto 606-8502, Japan

###### Abstract

Recent experiments on the Iron-based superconductor FeSe in a high magnetic
field have suggested the presence of both the fluctuation-induced vortex
liquid regime and a Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) vortex lattice. To
get a general picture on the magnetic phase diagram in type II superconductors
with strong superconducting (SC) fluctuation and strong paramagnetic pair-
breaking (PPB) such as FeSe, the vortex lattice melting curve $H_{m}(T)$ is
theoretically investigated in the situations where a FFLO state is expected to
occur. In general, PPB tends to narrow the vortex liquid regime intervening
between $H_{c2}(T)$ and $H_{m}(T)$. In particular, the vortex liquid regime is
found to rapidly shrink upon entering, by cooling, the temperature range in
which the FFLO state with a periodic modulation parallel to the magnetic field
is stable in the mean field theory. Based on the present results, the high
field SC phase diagrams of FeSe in the parallel and perpendicular field
configurations are discussed.

##  I Introduction

In type II superconductors, two kinds of spatially modulations are created in
the superconducting (SC) states by an applied magnetic field. One is the
vortex structure Abrikosov  protected by the flux quantization, which is a
topological condition in real space, and the other is a kind of Fulde-Ferrell-
Larkin-Ovchinnikov (FFLO) spatial modulation FF ; LO  which is usually
supported by a field-induced splitting between the up-spin and down-spin
portions of the Fermi surface. The former is a consequence of the orbital-pair
breaking effect of the magnetic field, while the latter is that of the Pauli
paramagnetic pair-breaking (PPB).

The presence of the vortices qualitatively changes the nature of the SC
fluctuation: The fluctuation changes the SC transition line (or, the upper
critical field) $H_{c2}(T)$ in the mean field theory to a crossover line
between the regions with strongly interacting SC fluctuations and with weakly
interacting fluctuations IOT89 ; FFH . The former region below the $H_{c2}(T)$
curve is often called as the vortex liquid regime. In clean limit, the true SC
transition line in a magnetic field is the position of the vortex lattice
melting in the field v.s. temperature ($H$-$T$) phase diagram and separates
the vortex liquid regime from the vortex lattice or solid FFH ; BG . So far,
any superconducting material, including the high $T_{c}$ cuprates, with a
broad vortex liquid regime has not shown PPB-induced SC phenomena. On the
other hand, most of searches for a FFLO state have been so far performed in
situations where the effects of fluctuating vortices are invisible Bianchi ;
Kenzel .

Recent experiments on the Iron-based quasi two-dimensional (Q2D)
superconductor FeSe have shown that, in both cases with a field parallel to
the SC layers (${\bf H}\perp c$) Kasa20  and a field perpendicular to the
layers (${\bf H}\parallel c$) Kasa21 , the fluctuation-induced vortex liquid
regime coexists with a PPB-induced novel high field SC phase in the same phase
diagram. To the best of our knowledge, this is the first material with a
possible FFLO vortex phase existing just below a well-defined vortex liquid
regime. One remarkable feature on the high field phase diagram of FeSe is that
the mean field $H_{c2}$-transition in FeSe is apparently continuous in both
field configurations Kasa20 ; Kasa21 ; Hardy  in contrast to the corresponding
ones in CeCoIn5 discussed repeatedly previously Bianchi ; Kenzel ; Kenzelmann
. On the other hand, the nature of the transition between the high field SC
phase and the low field vortex solid differ between the two field
configurations Kasa20 ; Kasa21 .

In the present work, we investigate how the position of the vortex lattice
melting of a Q2D type II superconductor is affected by strong PPB. Throughout
the present work, we focus on the configuration with an applied magnetic field
perpendicular to the basal plane of the superconductor AI03 , and the vortex
lattice structure is assumed to be the conventional hexagonal one. The melting
line will be examined in the following two manners. First, the elastic free
energy of the vortex lattice is derived by taking account of the possibility
of a formation of a FFLO spatial modulation parallel to the applied magnetic
field, and the Lindemann criterion on the melting line BG ; Moore  is derived
based on the obtained elastic energy. Second, the melting curve will also be
examined within the Ginzburg-Landau (GL) fluctuation analysis Hikami  by
comparing the free energy obtained by approaching from the normal state with
the free energy in the fluctuating vortex solid. In both of the two methods,
the obtained vortex liquid region becomes narrower with increasing PPB and, in
particular, upon entering the low temperature range in which the FFLO state is
stable in the mean field approximation. Qualitatively, this result is
consistent with the previous argument AI03  based on the Lindemann criterion
that the melting line should merge with the $H_{c2}$-line in the temperature
range where the mean field $H_{c2}$-transition is of first order. On the other
hand, it will also be clarified that the amplitude fluctuation makes the
validity of this argument vague.

The present paper is organized as follows. The vortex lattice melting field
$H_{m}(T)$ is derived according to the Lindemann criterion in sec.II and based
on comparison between the obtained free energies in sec.III. In sec.IV,
numerically obtained phase diagrams resulting from the methods in sec.II and
III are compared with each other. In sec.V, our results are summarized, and
their relevances to experimental phase diagrams of FeSe are discussed in
details. In Appendix, the microscopic derivation of the GL model is reviewed.

##  II Lindemann Criterion

The Ginzburg-Landau (GL) Hamiltonian we use in the main text of this paper
takes the form

$\displaystyle{\cal H}_{\rm GL}$ $\displaystyle=$ $\displaystyle N(0)\int
d^{3}{\bf
r}\biggl{[}\Delta^{*}(a_{0}+B(-\partial_{z}^{2})+C\partial_{z}^{4})\Delta+\frac{V_{4}}{2}|\Delta|^{4}$
(1) $\displaystyle+$ $\displaystyle\frac{V_{6}}{3}|\Delta|^{6}\biggr{]},$

where $N(0)$ is the electronic density of states per spin in the normal state.
In eq.(1), it was assumed that the order parameter $\Delta$ is already in its
lowest Landau level ($n=0$ LL) subspace, and that the mean field $H_{c2}(T)$
line is given by $a_{0}=0$ as far as $V_{4}$ is not negative. A microscopic
derivation of eq.(1) in the presence of PPB was performed previously AI03 ;
RI07 , and the details of the coefficients appearing in eq.(1) are given in
Appendix. In low enough fields and at higher temperatures where PPB is
negligible, the coefficients $B$ and $V_{4}$ are positive, while $C$ and
$V_{6}$ are negative and are conventionally assumed to be zero in such cases.
On the other hand, in the situations with strong PPB, the coefficients $C$ and
$V_{6}$ change their sign with decreasing the temperature and with increasing
the field, and, upon cooling further, $V_{4}$ and $B$ become negative in
higher fields. A negative $B$ implies that a FFLO state with a modulation
parallel to the field tends to form upon cooling. On the other hand, a
negative $V_{4}$ implies that the mean field $H_{c2}$ transition is of first
order, although, in real systems with SC fluctuation, this first order
transition never occurs in reality and is reflected just as a crossover.
Below, it will be examined how these PPB effects affect the vortex lattice
melting line $H_{m}(T)$.

The mean field vortex lattice solution is obtained in a conventional manner
Abrikosov . First, as far as any SC fluctuation is absent, the non-Gaussian
terms of eq.(1) can be rewritten in the manner

$\biggl{\langle}\frac{V_{4}}{2}|\Delta|^{4}+\frac{V_{6}}{3}|\Delta|^{6}\biggr{\rangle}_{s}=\frac{{\tilde{V}}_{4}}{2}(\langle|\Delta|^{2}\rangle_{s})^{2}+\frac{{\tilde{V}}_{6}}{3}(\langle|\Delta|^{2}\rangle_{s})^{3}$
(2)

where $\langle\,\,\,\,\rangle_{s}$ denotes the spatial average, $S$ is the
system area in the plane perpendicular to the field,
${\tilde{V}}_{4}=\beta_{\rm A}V_{4}$, and ${\tilde{V}}_{6}=\gamma_{\rm
A}V_{6}$ with

$\displaystyle\beta_{\rm A}$ $\displaystyle=$
$\displaystyle\frac{\langle|\Delta|^{4}\rangle_{s}}{(\langle|\Delta|^{2}\rangle_{s})^{2}}=1.1596,$
$\displaystyle\gamma_{\rm A}$ $\displaystyle=$
$\displaystyle\frac{\langle|\Delta|^{6}\rangle_{s}}{(\langle|\Delta|^{2}\rangle_{s})^{3}}=1.4230.$
(3)

Hereafter, as a possible PPB-induced spatial modulation of $\Delta$ developing
along the applied field direction, the helical phase modulation will also be
included. Then, the mean field solution with the in-plane triangular vortex
lattice structure and the out-of-plane helical phase modulation is given by
$\Delta_{0}=\alpha_{0}\varphi({\bf r}_{\perp}|0)e^{iq_{m}z}$ with
$q_{m}^{2}=\theta(-B)|B|/2C$, where $\varphi({\bf r}_{\perp}|0)$ is the
Abrikosov solution Abrikosov ; Moore ; Eilenberger ; RI591740

$\varphi({\bf
r}_{\perp}|0)=\sqrt{\frac{kr_{H}}{\pi^{1/2}}}\sum_{n=-\infty}^{\infty}\exp\biggl{[}-\frac{y^{2}}{2r_{H}^{2}}+{\rm
i}kn\biggl{(}x+\frac{\pi}{2k}n-iy\biggr{)}\biggr{]}$ (4)

constructed in the $n=0$ LL and in the Landau gauge ${\bf A}=-Hy{\hat{x}}$ and
satisfying the normalization condition $\langle|\varphi({\bf
r}|0)|^{2}\rangle_{s}=1$, $k=\pi^{1/2}3^{1/4}r_{H}^{-1}$,
$r_{H}=\sqrt{\phi_{0}/(2\pi H)}$, $\phi_{0}=\pi\hbar/|e|$ is the flux quantum,
and ${\bf r}_{\perp}=$ ($x$, $y$) denotes the 2D coordinate. The value of
$\alpha_{0}^{2}$ minimizing the free energy, eqs.(1) and (2), is given by

$\alpha_{0}^{2}=\frac{|{\tilde{V}}_{4}|}{2{\tilde{V}}_{6}}\biggl{(}-s_{4}+\sqrt{1-{\overline{a}}^{\prime}_{0}}\biggr{)},$
(5)

where

${\overline{a}}^{\prime}_{0}=\frac{4{\tilde{V}}_{6}}{{\tilde{V}}_{4}^{2}}a^{\prime}_{0}=\frac{4{\tilde{V}}_{6}}{{\tilde{V}}_{4}^{2}}\,(a_{0}-Cq_{m}^{2}),$
(6)

and $s_{4}=V_{4}/|V_{4}|$. The mean field SC transition at $H_{c2}$ is of
second order when $s_{4}>0$, while it is of first order when $s_{4}<0$. The
resulting $H_{c2}(T)$-line is given by $a^{\prime}_{0}=0$ for the former,
while it is given by ${\overline{a}}^{\prime}_{0}=3/4$ for the latter. Then,
the free energy density of the mean field solution becomes

$f_{\rm
MF}=-N(0)\frac{|{\tilde{V}}_{4}|^{3}}{12{\tilde{V}}^{2}_{6}}\biggl{[}\frac{1}{2}s_{4}{\overline{a}}^{\prime}_{0}+(1-{\overline{a}}^{\prime}_{0})^{3/2}-s_{4}(1-{\overline{a}}^{\prime}_{0})\biggr{]}.$
(7)

Next, the elastic energy of the vortex lattice will be considered. As far as
we restrict ourselves to the type II limit with no gauge field fluctuation
incorporated, the elastic energy of the vortex lattice is obtained as the
energy of the massless harmonic fluctuation within the $n=0$ LL around the
vortex lattice solution (4) Moore ; Eilenberger ; RI591740 . Its derivation in
the presence of PPB is sketched in Appendix. The resulting Hamiltonian of the
massless mode becomes

$\delta{\cal H}_{\rm ph}=\frac{1}{2}\sum_{q,{\bf
k}_{\perp}}(\rho_{s}q^{2}+\sigma_{s}q^{4}+C_{66}{\overline{k}}_{\perp}^{4})|\delta\chi(q,{\bf
k}_{\perp})|^{2},$ (8)

where $q$ is the wave number measured from $q_{m}$ in the $z$-direction AI03 ,

$\displaystyle\rho_{s}$ $\displaystyle=$ $\displaystyle
2N(0)\,\alpha_{0}^{2}\,|B|\,(1+\theta(-B)),$ $\displaystyle\sigma_{s}$
$\displaystyle=$ $\displaystyle 2N(0)\,C\,\alpha_{0}^{2}\,\theta(C),$ (9)

and

$C_{66}=2N(0)\,\alpha_{0}^{4}\,(0.119\,V_{4}+0.276\,\alpha_{0}^{2}\,V_{6}).$
(10)

The first term of eq.(10) coincides with the result in the previous work Moore
. In fact, in low enough fields ($H\ll H_{c2}(0)$) where PPB is negligible so
that $B$ approaches $\xi_{0}^{2}$, we have $\rho_{s}r_{H}^{-4}\simeq
H^{2}/(4\pi\lambda^{2}(T))$, and eq.(8) becomes the energy of the shear
elastic fluctuation in type II limit by assuming, as mentioned earlier,
$\sigma_{s}$ and $V_{6}$ to be zero:

$\delta{\cal H}_{\rm ph}\simeq\frac{1}{2}\sum_{q,{\bf
k}_{\perp}}\biggl{(}\frac{H^{2}}{4\pi\lambda^{2}(T)k_{\perp}^{2}}q^{2}+C_{66}k_{\perp}^{2}\biggr{)}|s_{\rm
T}(q,{\bf k}_{\perp})|^{2}$ (11)

where $r_{H}^{2}(\nabla\times{\hat{z}})\delta\chi$ was identified with the
transverse component ${\bf s}_{\rm T}$ of the vortex displacement field Moore
; RI591740 . Here, $\xi_{0}$ and $\lambda(0)$ are, respectively, the zero
temperature coherence length and penetration depth defined within the GL
theory in low fields.

Based on eq.(10), one might be afraid of whether $C_{66}$ approaches zero upon
cooling in the region where $V_{4}$ is negative. However, one can check that
$C_{66}$ remains positive even upon approaching the first order
$H_{c2}$-transition line from below.

The Lindemann criterion for determining the vortex lattice melting line
$H_{m}(T)$ becomes

$\langle s_{\rm T}^{2}\rangle=\int_{q}\int_{{\bf
k}_{\perp}}\frac{Tk_{\perp}^{2}r_{H}^{4}}{\rho_{s}q^{2}+\sigma_{s}q^{4}+C_{66}k_{\perp}^{4}r_{H}^{4}}=c_{\rm
L}^{2}r_{H}^{2},$ (12)

where the ${\bf k}_{\perp}$-integral is performed by setting the first
Brillouin Zone of the vortex lattice to be, for simplicity, circular. The
constant parameter $c_{\rm L}$ needs to be determined phenomenologically or
empirically. The $q$-integral can be performed analytically, and the Lindemann
criterion in the presence of PPB becomes

$\frac{T}{\sqrt{\rho_{s}C_{66}}}=c_{\rm
L}^{2}\frac{\phi_{0}}{H}\biggl{[}1+\sqrt{1+4\sqrt{\frac{C_{66}\sigma_{s}}{\rho_{s}^{2}}}}\biggr{]}.$
(13)

Based on the above-mentioned fact that, when the $H_{c2}$-transition is of
first order, $C_{66}$ does not vanish on approaching the $H_{c2}$-line from
below, it is clear that the equality in eq.(13) is not satisfied at low enough
$T$. It inevitably leads to the argument that, in the temperature range where
the mean field SC transition at $H_{c2}$ is of first order, the melting curve
$H_{m}(T)$ tends to merge with the $H_{c2}$-curve at a finite temperature, and
consequently that the vortex liquid regime in a superconducting material with
strong PPB tends to disappear at the temperature on cooling AI03 ; Hardy . The
validity of this argument based on the elastic theory will be discussed again
in sec.IV.

##  III Approach based on free energy evaluation

It is useful to compare the result on the melting line $H_{m}(T)$ in sec.II
with that following from a different approach in order to see to what extent
the Lindemann criterion in sec. II is reliable. In this section, we try to
obtain $H_{m}(T)$ based directly on calculating the fluctuation free energy
Hikami .

Imagine that one starts from the normal phase. The fluctuation propagator
${\cal D}_{p,q}=\langle|\Delta_{0,p,q}|^{2}\rangle$ of the $n=0$ LL modes
$\Delta_{0,p,q}$ in the Gaussian approximation is defined as

${\cal D}_{p,q}=\frac{T}{N(0)(a_{0}+Bq^{2}+Cq^{4})}.$ (14)

in the case of the GL Hamiltonian in eq.(1). Here, as in sec.II, the Landau
gauge was chosen for the vector potential. Due to the degeneracy in each LL,
the r.h.s. of eq.(14) is independent of the wavenumber $p$ defined in the
plane perpendicular to the field.

To extend this fluctuation propagator to the case with the mode-coupling
terms, the renormalized mass $\mu$ will be introduced as

$\mu=a_{0}+\Sigma.$ (15)

Then, eq.(14) is replaced by

${\cal D}_{p,q}=\frac{T}{N(0)(\mu+Bq^{2}+Cq^{4})}.$ (16)

In general, the self energy $\Sigma$ is constructed based on the perturbation
expansion with respect to the mode-coupling terms, and a non-Gaussian or
renormalized theory of the SC fluctuation is formulated IOT89 ; Hikami . In
contrast to the case Hikami ; Ruggeri  studied thoroughly so far with no PPB,
however, it is extremely cumbersome to study the large order behavior of the
perturbation series in the present case with PPB where the sixth order term
and the quartic gradient term need to be taken into account in the starting GL
free energy. Fortunately, the feature deep in the vortex liquid regime that
the amplitude $|\Delta|$ of the SC order parameter is well defined holds in
the present case with strong PPB. This feature is well approximated based on
the self-consistent Hartree-Fock approximation IOT90 ; AI17  in which $\mu$ is
independent of $q$. Below, the Hartree-Fock approximation will be used by
expecting the parameter dependences of the resulting phase diagram to be
comparable with those of the result of the Lindemann criterion. Then, the
renormalized mass $\mu$ is determined by the self-consistent equation (15)
with

$\displaystyle\Sigma$ $\displaystyle=$ $\displaystyle
V_{4}\frac{H}{N(0)\phi_{0}}\frac{T}{\sqrt{\mu(B+2\sqrt{C\mu})}}+\frac{3V_{6}}{2}\biggl{(}\frac{H}{N(0)\phi_{0}}\biggr{)}^{2}$
(17) $\displaystyle\times$ $\displaystyle\frac{T^{2}}{\mu(B+2\sqrt{C\mu})}.$

We note that, in the case with a negative $B$, i.e., in the temperature range
where the ordered phase in lower fields is a FF vortex lattice, not $\mu$ but

$\mu^{\prime}\equiv\mu-\frac{B^{2}}{4C}$ (18)

plays the role of the renormalized mass, reflecting the fact that the
corresponding bare mass is not $a_{0}$ but $a^{\prime}_{0}$ defined in eq.(6).

Next, to write down the expression of the fluctuation free energy density
$f_{>}$ in the normal phase, the exact expression on the mean squared average
of the SC order parameter $\langle|\Delta|^{2}\rangle$

$N(0)\langle|\Delta|^{2}\rangle=\frac{\partial f_{>}(a_{0})}{\partial a_{0}}$
(19)

will be used. Then, in the Hartree-Fock approximation, the free energy density
in $B>0$ case becomes

$\displaystyle f_{>}$ $\displaystyle=$
$\displaystyle\int_{\infty}^{a_{0}}d\varepsilon\frac{\partial
f_{>}(\varepsilon)}{\partial\varepsilon}+f_{0}$ (20) $\displaystyle=$
$\displaystyle
N(0)\int_{0}^{\mu}d\mu\,\langle|\Delta|^{2}\rangle-\int_{\infty}^{\mu}d\mu\frac{\partial\Sigma}{\partial\mu}\frac{\partial
f_{>}(a_{0})}{\partial a_{0}}$ $\displaystyle=$
$\displaystyle\frac{HT}{\phi_{0}}\biggl{[}\frac{\sqrt{B+2\sqrt{C\mu}}-\sqrt{|B|}}{\sqrt{C}}-\frac{V_{4}}{4}\biggl{(}\frac{HT}{N(0)\phi_{0}}\biggr{)}$
$\displaystyle\times$
$\displaystyle\frac{1}{\mu(B+2\sqrt{C\mu})}-\frac{V_{6}}{2}\biggl{(}\frac{HT}{N(0)\phi_{0}}\biggr{)}^{2}$
$\displaystyle\times$
$\displaystyle\frac{1}{[\mu(B+2\sqrt{C\mu})]^{3/2}}\biggr{]}.$

In the above analysis, the expression

$N(0)\int_{0}^{\infty}d\mu\langle|\Delta|^{2}\rangle$ (21)

has been chosen as the constant $f_{0}$ in the first line of eq.(20). Then, in
$C\to 0$ limit, eq.(20) is precisely the same as the corresponding expression
in Ref.Hikami .

The corresponding expression to eq.(20) in $B<0$ case is obtained in almost
the same manner, although
$N(0)\int_{0}^{\infty}d\mu^{\prime}\langle|\Delta|^{2}\rangle$ is chosen as
$f_{0}$ of the first line in eq.(20). Then, the free energy density in $B<0$
case is given simply by adding $HT(|B|/C)^{1/2}/\phi_{0}$ to eq.(20), that is,

$\displaystyle f_{>}$ $\displaystyle=$
$\displaystyle\frac{HT}{\phi_{0}}\biggl{[}\sqrt{\frac{4\mu^{\prime}}{2\sqrt{C\mu}+|B|}}-\frac{V_{4}}{16}\biggl{(}\frac{HT}{N(0)\phi_{0}}\biggr{)}\frac{|B|+2\sqrt{C\mu}}{C\mu\mu^{\prime}}$
(22) $\displaystyle-$
$\displaystyle\frac{V_{6}}{16}\biggl{(}\frac{HT}{N(0)\phi_{0}}\biggr{)}^{2}\biggl{[}\frac{|B|+2\sqrt{C\mu}}{C\mu\mu^{\prime}}\biggr{]}^{3/2}\,\biggr{]}.$

In turn, the free energy density of the vortex lattice state will be examined
Hikami . To do this, the free energy arising from the harmonic excitations
around the mean field vortex lattice will be added to $f_{\rm MF}$ given in
eq.(7). However, it is a well established fact within the approach based on
the critical SC fluctuation that the shear elastic mode of the vortex lattice,
i.e., the Goldstone mode, is a smaller correction compared with that of the
amplitude fluctuation. Therefore, the free energy term resulting from eq.(8)
can be neglected below. The Hamiltonian $\delta{\cal H}_{\rm amp}$ on the
Gaussian amplitude fluctuation is presented in Appendix. In the same way as
the derivation of the first term of $f_{>}$ in sec.II, the fluctuation
correction $\delta f$ to $f_{\rm MF}$ resulting from $\delta{\cal H}_{\rm
amp}$ is obtained in the form

$\delta
f=\frac{HT}{\sqrt{2C}\phi_{0}}\biggl{(}\,\sqrt{B+2\sqrt{r_{<}C}}-\sqrt{B}\,\biggr{)},$
(23)

for $B>0$, and

$\delta f=\frac{HT}{\sqrt{2C}\phi_{0}}\sqrt{2|B|+2\sqrt{r_{<}C}}$ (24)

for $B<0$, where

$r_{<}\simeq\frac{{\tilde{V}}_{4}^{2}}{{\tilde{V}}_{6}}\biggl{(}1-{\overline{a}}^{\prime}_{0}-s_{4}\sqrt{1-{\overline{a}}^{\prime}_{0}}\biggr{)}.$
(25)

In this manner, the free energy density based on the approach from lower
temperatures is given by $f_{<}=f_{\rm MF}+\delta f$.

After all, the melting line $H_{m}(T)$ is determined by the equality
$f_{>}=f_{<}$ within the approach in this section.

##  IV Results on Magnetic Phase Diagram

Figure 1: (Color online) Field v.s. temperature ($H$-$T$) phase diagrams
obtained (a) from eq.(13) and (b) from the approach in sec.III. The used
parameter values are $\alpha_{\rm M}=0.25$ and ${\tilde{N}}_{0}=0.1$. The red
dotted curve with red square symbols denotes the crossover line corresponding
to the mean field second order SC transition line $H_{c2}(T)$, while the blue
solid curve denotes the vortex lattice melting line which is the genuine SC
transition line in clean limit.

In this section, the melting transition lines of the $n=0$ LL vortex lattice
following from the two methods explained in the preceding sections will be
compared with each other to obtain a generic picture on the magnetic ($H$-$T$)
phase diagram of the type II superconductors with moderately strong PPB. To
make understanding of the parameter dependences of the phase diagram easier,
it will be useful to rewrite the formula eq.(13) based on the Lindemann
criterion in terms of dimensionless variables in the form

$\frac{T}{4\pi
T_{c0}}\frac{\xi_{0}^{2}}{r_{H}^{2}}={\tilde{N}}_{0}\biggl{(}\frac{c_{\rm
L}\alpha_{0}}{T_{c0}}\biggr{)}^{2}g\,\biggl{[}|b|^{1/2}+\sqrt{|b|+4|c|^{1/2}g}\biggr{]},$
(26)

where

$g=\sqrt{0.119\alpha_{0}^{2}{\tilde{V}}_{4}+0.276\alpha_{0}^{4}{\tilde{V}}_{6}},$
(27)

and

${\tilde{N}}_{0}=N(0)T_{c0}\xi_{0}^{2}\xi_{0,\parallel}=\frac{0.003}{\sqrt{\rm
Gi}}$ (28)

expressed in terms of the Ginzburg number

$\displaystyle{\rm Gi}$ $\displaystyle=$ $\displaystyle
2\biggl{(}\frac{7\zeta(3)}{64\pi^{3}N(0)T_{c0}\xi_{0}^{2}\xi_{0,\parallel}}\biggr{)}^{2}$
(29) $\displaystyle=$ $\displaystyle
2\biggl{(}\frac{7\zeta(3)}{32\pi^{2}}\,\frac{(\lambda(0))^{2}}{\Lambda(T_{c0})\xi_{0,\parallel}}\biggr{)}^{2}\,\frac{|\Delta(0)|^{4}}{T_{c0}^{4}}$

measuring the strength of the thermal fluctuation is the scaled DOS, $\xi_{0}$
and $\xi_{0,\parallel}$ are the in-plane and out-of-plane coherence lengths of
a Q2D material, $|\Delta(0)|$ is the zero temperature energy gap,
$\Lambda(T)=\phi_{0}^{2}/(16\pi^{2}T)$ is the thermal length FFH , and the
dimensionless coefficients, $b=B/\xi_{0,\parallel}^{2}$ and
$c=C/\xi_{0,\parallel}^{4}$, are given in (37) in Appendix. Based on eq.(26),
$H_{m}(T)$ will be discussed hereafter as a function of the fluctuation
strength $1/{\tilde{N}}_{0}$ and the PPB strength, i.e., the Maki parameter
Maki

$\alpha_{\rm M}=\frac{\mu H_{c2}^{({\rm orb})}(T=0)}{2\pi T_{c0}}$ (30)

which is incorporated in the dimensionless GL coefficients $a_{0}$, $b$, $c$,
$\alpha_{0}/T_{c0}$, $T_{c0}^{2}V_{4}$, and $T_{c0}^{4}V_{6}$, where $\mu H$
is the Zeeman energy for a single quasiparticle, and $H_{c2}^{({\rm orb})}(T)$
is the $H_{c2}(T)$-line in ${\bf H}\parallel c$ in the absence of PPB.
Similarly, the contributions to the free energy density introduced in sec.III,
$f_{\rm MF}$, $f_{>}$, and $\delta f$, are also described in terms of the
dimensionless GL coefficients parameterized by ${\tilde{N}}_{0}$ and
$\alpha_{\rm M}$. In all of the $H$-$T$ phase diagrams to be discussed below,
the temperature $T$ and the field strength $H$ are expressed in the units of
$T_{c0}$ and the Pauli-limiting field $T_{c0}/\mu=H_{c2}^{({\rm
orb})}(T=0)/(2\pi\alpha_{\rm M})$, respectively.

Hereafter, the parameter values $\alpha_{\rm M}=0.25$ and $0.75$ will be used
together with ${\tilde{N}}_{0}=0.1$ and $1.0$. Corresponding to the Gi-value
of FeSe noted in Ref.Koshelev , we have

${\tilde{N}}_{0}=0.12.$ (31)

Figure 2: (Color online) Field v.s. temperature ($H$-$T$) phase diagrams
obtained (a) from eq.(13) and (b) from the approach in sec.III. The used
parameter values are $\alpha_{\rm M}=0.75$ and ${\tilde{N}}_{0}=1.0$. The
$H_{c2}(T)$-curve (red square symbols) consists of the mean field second order
transition line (red dotted curve) at which ${\overline{a}}^{\prime}_{0}=0$
above $T^{*}$ (i.e., when $V_{4}>0$) and the mean field first order transition
line (the black-dotted one) at which ${\overline{a}}^{\prime}_{0}=3/4$ below
$T^{*}$. The green solid curve at which the coefficient $B$ changes its sign
is the transition line within $n=0$ LL between the conventional Abrikosov
vortex lattice and the FFLO vortex lattice at lower temperatures. The melting
line is expressed by the blue solid curve.

Figures 1 and 2 are different examples of comparison between the two melting
lines obtained from the methods introduced in the preceding two sections under
a fixed set of ${\tilde{N}}_{0}$ and $\alpha_{\rm M}$ values. In both of
Figs.1 and 2, the melting curves obtained in terms of the two methods nearly
coincide with each other when using the Lindemann constant $c_{\rm L}=0.414$,
suggesting that the parameter dependences determining $H_{m}(T)$ are similar
to each other between the two methods. Hereafter, the value $c_{\rm L}=0.414$
will be commonly used in obtaining results based on the Lindemann criterion.
It will be seen later that such agreement is also seen in the case with
stronger PPB and stronger fluctuation at least outside the FFLO temperature
range in which the FFLO state occurs in the mean field theory (see Fig.3
below).

In Fig.1, a strong fluctuation strength ${\tilde{N}}_{0}^{-1}=10$ and a
relatively weaker PPB strength $\alpha_{\rm M}=0.25$ are used, and, as in the
case with no PPB, $H_{m}(T)$ is concave, i. e., a curve with a positive
curvature in the $H$-$T$ phase diagram and approximately obeys the $n=0$ LL
scaling IOT89 ; Moore ; Hikami  at least in lower fields where $B>0$ and
$C=0$. At low enough temperatures where $B$ ($>0$) is small, and $C>0$, the
$n=0$ LL scaling is not satisfied any longer. Nevertheless, the deviation from
the $n=0$ LL scaling seems to be unexpectedly small. In contrast, in Fig.2
with ${\tilde{N}}_{0}^{-1}=1.0$ and $\alpha_{\rm M}=0.75$, the vortex liquid
regime is narrow, and $H_{m}(T)$ begins to follow the functional form of the
$H_{c2}(T)$ curve. That is, $H_{m}(T)$ is convex and has a negative curvature,
except in the close vicinity of $T_{c0}$, in the $H$-$T$ phase diagram.

Note that the $H_{m}(T)$-curve in Fig.2 apparently merges with the
$H_{c2}(T)$-curve at a higher temperature than the FFLO temperature range.
This tendency implying a shrinkage of the vortex liquid regime at a finite
temperature is in agreement with the argument noted at the end of sec.II AI03
. As is seen in Fig.3, however, this tendency becomes unclear when the
fluctuation is stronger.

Figure 3: (Color online) Field v.s. temperature ($H$-$T$) phase diagrams
obtained (a) from eq.(13) and (b) from the approach in sec.III. The content of
the figures are the same as that of Fig.2 (a) and (b) except the use of the
parameter values $\alpha_{\rm M}=0.75$ and ${\tilde{N}}_{0}=0.1$ in these
figures. The portion (green dotted curve) of the mean field FFLO transition
line in the vortex liquid regime is just a crossover line.

To understand the ${\tilde{N}}_{0}$ and $\alpha_{\rm M}$ dependences of the
phase diagram in more details, the results on the phase diagram following from
the two methods are shown in Fig.3, where the values ${\tilde{N}}_{0}^{-1}=10$
and $\alpha_{\rm M}=0.75$ are used. These two figures show typical phase
diagrams in the case where both the fluctuation and PPB are moderately strong.
One main feature is that, in spite of the strong fluctuation, the
$H_{m}(T)$-curve is convex over most of the field range in the $H$-$T$ phase
diagram. On the other hand, concave portions of $H_{m}(T)$ are seen close to
$T_{c0}$ and in the FFLO temperature range. The former occurring in low enough
fields where PPB weakly contributes is a consequence of the $n=0$ LL scaling
IOT89 ; Moore ; Hikami  $T-T_{c0}\sim(TH)^{2/3}$, while the latter reflects
the shrinkage of the vortex liquid due to PPB AI03 .

By comparing Figs.2 and 3 with each other, it is easily found that, at a fixed
$\alpha_{\rm M}$, an increase of ${\tilde{N}}_{0}$ leads to a shrinkage of the
vortex liquid regime while the convex $H_{m}(T)$-curve is kept. This feature
suggests that, in superconductors with weak fluctuation, it is not easy to
distinguish $H_{m}(T)$ from $H_{c2}(T)$ through experimental data.

On the other hand, by comparing Figs.1 and 3 with each other, it is found that
an increase of PPB makes the vortex liquid regime narrower. In particular, as
mentioned above, it is commonly seen that the $H_{m}(T)$ curve becomes convex,
like $H_{c2}(T)$, in the $H$-$T$ phase diagram, reflecting an enhanced role of
PPB. Since ${\tilde{N}}_{0}^{-1}$ measures the strength of the fluctuation in
zero field, the feature mentioned above implies that, in systems with
moderately strong PPB, the actual fluctuation strength tends to be
underestimated through experimental data in finite fields by, for instances,
identifying the irreversibility line on which the resistivity vanishes with
the $H_{c2}(T)$-line, because the melting curve lies quite close to the
irreversibility line in most cases.

Next, let us discuss the fate of the melting line in the FFLO temperature
range. In Fig.3 (a) obtained based on the Lindemann criterion, the
$H_{m}(T)$-line suddenly begins to approach the $H_{c2}(T)$-line on entering
the FFLO temperature range. The origin of this sharp change seems to consist
in the change in the coefficient of the $q^{2}$ term upon entering the FFLO
temperature range by cooling (see eq.(9)). As seen in Fig.3 (b) obtained based
on the free energy approach explained in sec.III, however, this change in
$H_{m}(T)$ upon entering the FFLO range seems to become unclear as the
amplitude fluctuation is incorporated. Therefore, the fate of $H_{m}(T)$ in
the FFLO temperature range is sensitive to the details of its derivation, and
it is not sufficiently understood at present whether $H_{m}(T)$ truly merges
with $H_{c2}(T)$ in systems with strong enough fluctuation.

##  V Summary and Discussion

In the present work, the vortex lattice melting curve in the type II
superconductor with moderately strong paramagnetic pair-breaking (PPB) has
been theoretically examined by assuming the vortex lattice to have the
familiar hexagonal symmetry and hence to be described by the lowest ($n=0$) LL
modes of the SC order parameter. The present result extends the previous works
Moore ; Houghton  to the cases with PPB.

Below, let us discuss the magnetic phase diagrams of FeSe Kasa20 ; Kasa21 ;
Hardy  as an example of application of the results in the preceding sections.
The field configuration assumed in the present work corresponds to FeSe in
${\bf H}\parallel c$ where a nearly linear $H_{m}(T)$-line Hardy  and the
presence of a high field SC (HFSC) phase Kasa21  were found. The nearly linear
$H_{m}(T)$-curve is easily understood based on Fig.1 and Fig.3. Clearly, the
deviation from the concave melting curve suggesting the $n=0$ LL scaling Moore
; RI591740  is due to moderately strong PPB, and the region in which the $n=0$
LL scaling is visible is limited to the low field range close to $T_{c0}$. It
has been argued elsewhere AI19  that the low field behavior of the melting
line is also affected by the strong-coupling effect due to the SC fluctuation
itself in a system close to the so-called BCS-BEC crossover regime.

Figure 4: (Color online) Field dependences of the mean-squared amplitude
$\langle|\Delta|^{2}\rangle$ of the SC order parameter at the temperatures
$T/T_{c0}=0.2$ (blue), $0.15$ (green), $0.10$ (purple), $0.05$ (red), and
$0.01$ (black dotted curve) following from eqs.(15) and (17) in terms of
$\alpha_{\rm M}=0.8$ and ${\tilde{N}}_{0}=0.12$. Reflecting the fact that,
when $\alpha_{\rm M}=0.8$, the $H_{c2}(T)$-line is the first order transition
line in the mean field approximation in $t<0.35$, $\langle|\Delta|^{2}\rangle$
at $t=0.01$ where the thermal fluctuation is extremely weak shows a nearly
discontinuous change at $H_{c2}(0)$. Inclusion of the quantum fluctuation
neglected here would broaden even the $t=0.01$ curve.

On the other hand, extensive discussion is needed to understand the HFSC phase
in ${\bf H}\parallel c$ Kasa21 . By taking account of the fact that, in ${\bf
H}\parallel c$ where the fluctuation effect may become relatively strong, the
effective strength of PPB does not become extremely strong, we have at most
two candidates of such a high field SC phase induced by PPB. One is the second
lowest ($n=1$) LL vortex lattice Klein ; Allan ; Matsu . Within the analysis
based on the weak-coupling BCS approach for a single-band electronic model,
the low temperature phase induced by PPB in clean limit inevitably becomes the
$n=1$ LL vortex lattice (see Fig.6 (a) in Ref.AI03 ). However, the HFSC phase
in FeSe in ${\bf H}\parallel c$ case is not the $n=1$ LL state, judging from
the experimental fact Kasa21 ; Kasa14  that the field-induced structural
transition between the conventional vortex solid in lower fields and the HFSC
one is continuous. Then, the HFSC phase in ${\bf H}\parallel c$ should be
identified with one of the FFLO vortex solids, describable in the $n=0$ LL,
with a spatial modulation parallel to the applied field AI03 . According to
Fig.2 in Ref.AI15 , the $n=1$ LL vortex state may be destabilized even at zero
temperature in the cases of two-band electronic models with an $\alpha_{\rm
M}$-value of order unity so that it is possible that the HFSC phase in a two-
band system in clean limit is indeed a $n=0$ LL vortex lattice with a spatial
modulation parallel to ${\bf H}$. On the other hand, in the case of a $n=0$
vortex lattice described by the GL model, the type of the FFLO spatial
modulation is found to depend on the nature of the $H_{c2}$-transition at low
temperatures Agter ; Nakashima : When the mean field $H_{c2}$-transition is of
second order, the HFSC phase should be a hybrid of the triangular vortex
lattice and the phase-modulated FF state FF  and cannot become the triangular
vortex lattice with the LO-like periodic amplitude-modulation LO  parallel to
${\bf H}$. The triangular vortex lattice with the LO-like spatial modulation
parallel to ${\bf H}$ becomes the HFSC phase only when the mean field
$H_{c2}$-transition is of first order. Then, one might wonder if the fact
Kasa21  mentioned in sec.I that the $H_{c2}$-transition in FeSe in ${\bf
H}\parallel c$ seems to be continuous contradict the observation of a nodal
plane perpendicular to ${\bf H}$ at the sample surface strongly suggesting
that the HFSC phase in FeSe in ${\bf H}\parallel c$ should be the LO vortex
lattice Kasa21 . However, we find that, once the strong SC fluctuation in FeSe
is taken into account, there may be no such contradiction: By using the
relations eqs.(15), (17), (31), (37), and the data of $V_{4}$ and $V_{6}$ for
$\alpha_{\rm M}=0.8$ given in Fig.5 of Ref.AI03 , we obtain Fig.4 expressing
the field dependence of the mean-squared SC order parameter
$\langle|\Delta|^{2}\rangle$ at each temperature, $t=T/T_{c0}=0.2$, $0.15$,
$0.1$, $0.05$, and $0.01$. We note that the mean field $H_{c2}$-transition is
of first order, i.e., $V_{4}<0$, in $t\leq 0.35$ for the set of the parameters
used in Fig.4. Nevertheless, the ${\tilde{N}}_{0}$-value of eq.(31) makes the
discontinuous change of $\langle|\Delta|^{2}\rangle$ in the mean field theory
at $H_{c2}(0)$ broad enough in $t\geq 0.05$ (compare Fig.4 with Figs.9 and 11
in Ref.AI03 ). Further, the quantum SC fluctuation has not been taken into
account in the present analysis. Inclusion of the quantum fluctuation would
broaden even the $\langle|\Delta|^{2}\rangle$-curve at $t=0.01$ in Fig.4. By
comparing the curves in Fig.4 with the heat capacity data in Ref.Kasa21 , we
find that it is difficult to determine the nature of the mean field
$H_{c2}$-transition from the experimental data of real systems with strong
fluctuation in which the $H_{c2}$ is merely a continuous crossover line. Thus,
the present theory does not contradict the conclusion in Ref.Kasa21
identifying the HFSC phase in FeSe in ${\bf H}\parallel c$ as the LO vortex
lattice.

Finally, the magnetic phase diagram of FeSe in ${\bf H}\perp c$ will be
briefly discussed. Although the present work was performed by assuming the
configuration with a field perpendicular to the basal plane in a Q2D system,
the results in sec.IV should be qualitatively applicable even to the ${\bf
H}\perp c$ case. First, the $H_{m}(T)$-curve in Fig.3 is quite similar to the
irreversibility line in Ref.Kasa20  and to the melting line estimated in
Ref.Hardy  in that the $H_{m}(T)$-curve is convex in the $H$-$T$ phase diagram
in spite of showing a broad vortex liquid regime. However, we expect the FFLO
phase in the high $H$ and low $T$ corner in Fig.3 to, in this ${\bf H}\perp c$
case, have been replaced by the $n=1$ LL vortex lattice. In Ref.NNI , we have
shown that a peculiar field dependence of the resistive behavior at low enough
temperatures in ${\bf H}\perp c$ is qualitatively consistent with the
resistivity curve resulting from the quantum SC fluctuation in the $n=1$ LL in
the case with moderately strong PPB. It strongly suggests that the HFSC phase
in ${\bf H}\perp c$ should be the $n=1$ LL vortex solid Klein ; Allan ; Matsu
. Then, the vortex liquid controlled by the $n=1$ LL modes of the SC order
parameter should be present just above the melting line. Theoretical
description of such a novel vortex state should be left for a future work.

##  VI Acknowledgement

We thank Naratip Nunchot for his help on numerical analysis and for reading
the original manuscript and Yuji Matsuda and Shigeru Kasahara for discussions
on their experimental data. The present work was supported by KAKENHI
(No.21K03468) from Japan Society for the Promotion of Science.

##  VII Appendix

Our starting model for deriving the Ginzburg-Landau (GL) Hamiltonian is the
simplest BCS Hamiltonian with a single electronic band

$\displaystyle{\cal H}$ $\displaystyle=$ $\displaystyle\sum_{\sigma=\pm 1}\int
d^{3}r\varphi^{\dagger}_{\sigma}({\bf
r})\biggl{[}\frac{\hbar^{2}}{2m}\biggl{(}-{\rm
i}\nabla+\frac{\pi}{\phi_{0}}{\bf
A}\biggr{)}^{2}-I\sigma\biggr{]}\varphi_{\sigma}({\bf r})$ (32)
$\displaystyle-$ $\displaystyle|g|\sum_{\bf q}\Psi^{\dagger}({\bf q})\Psi({\bf
q}),$

where $\phi_{0}=\pi\hbar/|e|$ is the flux quantum, $I=\mu H$ is the Zeeman
energy, $|g|$ is the attractive interaction strength, and

$\Psi({\bf q})=\frac{1}{2}\sum_{\bf p}\sum_{\sigma=\pm 1}\,\sigma\,w_{\bf
p}\,c_{{-{\bf p}+{\bf q}/2},-\sigma}\,c_{{\bf p}+{\bf q}/2,\sigma}$ (33)

is the pair-field operator expressed by a spin-singlet pairing function
$w_{\bf p}$ and $c_{{\bf p},\sigma}$ which is the Fourier transform of the
electron operator $\varphi_{\sigma}({\bf r})$. For simplicity, the $s$-wave
paired case with $w_{\bf p}=1$ will be assumed. The following GL Hamiltonian
is obtained from the electronic model (32) through an extension to Q2D case:

${\cal H}_{\rm GL}=N(0)\int d^{3}{\bf
r}\biggl{[}\Delta^{*}(a_{0}+B(-\partial_{z}^{2})+C\partial_{z}^{4})\Delta+\frac{V_{4}}{2}|\Delta|^{4}+\frac{V_{6}}{3}|\Delta|^{6}\biggr{]},$
(34)

where the order parameter $\Delta$ is assumed to be in the $n=0$ LL subspace,
and the applied field is assumed to be perpendicular to the basal plane. The
coefficient $a_{0}$ is given by

$a_{0}={\rm ln}(t)+\int_{0}^{\infty}d\rho\biggl{(}\frac{2\pi t}{{\rm
sinh}(2\pi t\rho)}-f(\rho)\exp(-|\mu|^{2}\rho^{2}/2)\biggr{)},$ (35)

where $t=T/T_{c0}$, $T_{c0}$ is the zero field SC transition temperature, and

$f(\rho)=\frac{2\pi t}{{\rm sinh}(2\pi t\rho)}\,{\rm
cos}\biggl{(}2\frac{I}{T_{c0}}\rho\biggr{)}.$ (36)

The coefficients of the gradient terms are given by $B=b\xi_{0,\parallel}^{2}$
and $C=c\xi_{0,\parallel}^{4}$, where

$\displaystyle b$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\rho\rho^{2}f(\rho)\exp(-\rho^{2}|\mu|^{2}/2),$
$\displaystyle c$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\int_{0}^{\infty}d\rho\rho^{4}f(\rho)\exp(-\rho^{2}|\mu|^{2}/2),$
(37)

where $\xi_{0,\parallel}$ is the out-of-plane coherence length at zero
temperature, and $\mu=\sqrt{2}\pi\xi_{0}({\hat{p}}_{x}+{\rm
i}{\hat{p}}_{y})/r_{H}$.

The limitation to the $n=0$ LL modes imply that the mode coupling terms, the
fourth-order and six-order terms in eq.(34), are spatially nonlocal. As
indicated in Ref.AI03 , however, this nonlocality seems to be safely
negligible for most purposes, and, for simplicity, the local forms of the mode
coupling terms were assumed above. In the case of a layered system with a
cylindrical Fermi surface, the coefficients $V_{4}$ and $V_{6}$ have been
derived in Ref.AI03 ; RI07  and, in clean limit, are given by

$\displaystyle V_{4}$ $\displaystyle=$ $\displaystyle
3\,T_{c0}^{-2}\int\Pi_{i=1}^{3}d\rho_{i}f\biggl{(}\sum_{j=1}^{3}\rho_{j}\biggr{)}\biggl{\langle}\exp\biggl{(}-\frac{1}{2}\biggl{(}R_{14}-\frac{1}{2}R_{24}\biggr{)}\biggr{)}\,{\rm
cos}(I_{4})\biggr{\rangle}_{\rm FS},$ $\displaystyle V_{6}$ $\displaystyle=$
$\displaystyle-15T_{c0}^{-4}\int\Pi_{j=1}^{5}d\rho_{j}f\biggl{(}\sum_{j=1}^{5}\rho_{j}\biggr{)}\biggl{\langle}\exp\biggl{(}-\frac{1}{2}(R_{16}+R_{26})\biggr{)}\,{\rm
cos}(I_{6})\biggr{\rangle}_{\rm FS},$ (38)

where

$\displaystyle R_{14}$ $\displaystyle=$
$\displaystyle|\mu|^{2}(\sum_{j=1}{3}\rho_{j}^{2}+\rho_{2}(\rho_{3}+\rho_{1})),$
$\displaystyle R_{24}$ $\displaystyle=$ $\displaystyle{\rm
Re}(\mu^{2})(\rho_{2}^{2}+(\rho_{3}-\rho_{1})^{2})$ $\displaystyle I_{4}$
$\displaystyle=$ $\displaystyle\frac{{\rm
Im}(\mu^{2})}{4}(\rho_{2}^{2}-(\rho_{3}-\rho_{1})^{2})$ $\displaystyle R_{16}$
$\displaystyle=$
$\displaystyle|\mu|^{2}\biggl{(}e_{1}+e_{2}+e_{3}+\frac{2}{3}e_{4}e_{5}\biggr{)}$
$\displaystyle R_{26}$ $\displaystyle=$ $\displaystyle{\rm
Re}(\mu^{2})\biggl{(}e_{1}+e_{2}+e_{3}-\frac{e_{4}^{2}+e_{5}^{2}}{3}$
$\displaystyle-$ $\displaystyle\frac{2}{3}(e_{6}+e_{7}+e_{8}+e_{9})\biggr{)}$
$\displaystyle I_{6}$ $\displaystyle=$ $\displaystyle\frac{{\rm
Im}(\mu^{2})}{4}\biggl{(}e_{1}+e_{2}-e_{3}-\frac{e_{4}^{2}-e_{5}^{2}}{3}$ (39)
$\displaystyle-$ $\displaystyle\frac{2}{3}(e_{6}+e_{7}-e_{8}-e_{9})\biggr{)}$

$\displaystyle e_{1}$ $\displaystyle=$
$\displaystyle(\rho_{3}+\rho_{5})^{2}+(\rho_{3}+\rho_{4})^{2},$ $\displaystyle
e_{2}$ $\displaystyle=$ $\displaystyle(\rho_{1}+\rho_{4}+\rho_{5})^{2},$
$\displaystyle e_{3}$ $\displaystyle=$
$\displaystyle\rho_{3}^{2}+\rho_{4}^{2}+(\rho_{2}-\rho_{5})^{2},$
$\displaystyle e_{4}$ $\displaystyle=$
$\displaystyle\rho_{1}+2(\rho_{3}+\rho_{4}+\rho_{5}),$ $\displaystyle e_{5}$
$\displaystyle=$ $\displaystyle\rho_{2}-\rho_{3}-\rho_{4}-\rho_{5},$
$\displaystyle e_{6}$ $\displaystyle=$
$\displaystyle(\rho_{4}-\rho_{5})^{2}+(\rho_{1}+\rho_{5}-\rho_{3})^{2},$
$\displaystyle e_{7}$ $\displaystyle=$
$\displaystyle(\rho_{+}\rho_{4}-\rho_{3})^{2},$ $\displaystyle e_{8}$
$\displaystyle=$
$\displaystyle(\rho_{3}-\rho_{4})^{2}+(\rho_{2}+\rho_{3}-\rho_{5})^{2},$
$\displaystyle e_{9}$ $\displaystyle=$
$\displaystyle(\rho_{2}+\rho_{4}-\rho_{5})^{2}.$ (40)

To obtain the dispersion relations of the normal modes of the Gaussian
fluctuation around $\Delta_{0}({\bf r})$, the $q$-dependent terms, we follow
the Eilenberger’s analysis Eilenberger  to represent the total pair-field in
the form $\Delta=e^{{\rm i}q_{m}z}(\alpha_{0}\varphi({\bf
r}|0)+a_{+}\varphi({\bf r}|{\bf r}_{0})e^{{\rm i}qz}+a_{-}\varphi({\bf
r}|-{\bf r}_{0})e^{-{\rm i}qz})$. Then, in a situation with $b<0$, the terms
harmonic with respect to $a_{\pm}$ in ${\cal H}$ take the form

$\displaystyle\frac{\delta{\cal H}}{N(0)}$ $\displaystyle=$
$\displaystyle[a^{\prime}_{0}+Cq^{2}(q^{2}+4q_{m}^{2})+\alpha_{0}^{2}(2{\tilde{V}}_{4}{\overline{\xi}}_{d}({\bf
r}_{0})+3{\tilde{V}}_{6}\alpha_{0}^{2}{\overline{\eta}}_{d}({\bf
r}_{0}))](|a_{+}|^{2}+|a_{-}|^{2})+4Cq_{m}q^{3}(|a_{+}|^{2}-|a_{-}|^{2})$ (41)
$\displaystyle+$
$\displaystyle\alpha_{0}^{2}[({\tilde{V}}_{4}{\overline{\xi}}_{a}({\bf
r}_{0})+2{\tilde{V}}_{6}\alpha_{0}^{2}{\overline{\eta}}_{a}({\bf
r}_{0}))a_{+}a_{-}+{\rm c.c.}],$

where

$\displaystyle{\overline{\xi}}_{d}({\bf r}_{0})$ $\displaystyle=$
$\displaystyle\beta_{\rm A}^{-1}\langle|\varphi({\bf r}|0)\varphi({\bf r}|{\bf
r}_{0})|^{2}\rangle_{s},$ $\displaystyle{\overline{\xi}}_{a}({\bf r}_{0})$
$\displaystyle=$ $\displaystyle\beta_{\rm A}^{-1}\langle(\varphi^{*}({\bf
r}|0))^{2}(\varphi({\bf r}|{\bf r}_{0}))^{2}\rangle_{s},$
$\displaystyle{\overline{\eta}}_{d}({\bf r}_{0})$ $\displaystyle=$
$\displaystyle\gamma_{\rm A}^{-1}\langle|\varphi({\bf r}|0)\varphi({\bf
r}|{\bf r}_{0})|^{2}|\varphi({\bf r}|0)|^{2}\rangle_{s},$
$\displaystyle{\overline{\eta}}_{a}({\bf r}_{0})$ $\displaystyle=$
$\displaystyle\gamma_{\rm A}^{-1}\langle(\varphi^{*}({\bf
r}|0))^{2}(\varphi({\bf r}|{\bf r}_{0}))^{2}|\varphi({\bf
r}|0)|^{2}\rangle_{s}.$ (42)

By performing the ”Bogoliubov transformation”, the diagonalized form of
eq.(41) becomes $\delta{\cal
H}=N(0)(E_{+}|{\tilde{a}}_{+}|^{2}+E_{-}|{\tilde{a}}_{-}|^{2})$, where

$\displaystyle E_{\pm}$ $\displaystyle=$
$\displaystyle\frac{{\tilde{V}}_{4}^{2}}{4{\tilde{V}}_{6}}\biggl{(}-s_{4}+\sqrt{1-{\overline{a}}^{\prime}_{0}}\biggr{)}\biggl{[}4s_{4}{\overline{\xi}}_{d}+3\biggl{(}-s_{4}$
(43) $\displaystyle+$
$\displaystyle\sqrt{1-{\overline{a}}^{\prime}_{0}}\biggr{)}{\overline{\eta}}_{d}\pm
2\biggl{|}s_{4}{\overline{\xi}}_{a}+\biggl{(}-s_{4}+\sqrt{1-{\overline{a}}^{\prime}_{0}}\biggr{)}{\overline{\eta}}_{a}\biggr{|}\biggr{]}$
$\displaystyle+$ $\displaystyle a^{\prime}_{0}+Cq^{4}+2|B|q^{2}.$

In obtaining eq.(43), O($q^{6}$) terms were neglected. Here,
$\sqrt{2}a_{-}=\sigma e^{-{\rm i}\gamma/2}{\tilde{a}}^{*}_{\sigma}$,
$\sqrt{2}a_{+}=e^{-{\rm i}\gamma/2}{\tilde{a}}_{\sigma}$ ($\sigma=\pm$), and
$\gamma$ is the phase of
$s_{4}{\overline{\xi}}_{a}+(-s_{4}+\sqrt{1-{\overline{a}}^{\prime}_{0}}){\overline{\eta}}_{a}$.
The corresponding expression in the case with a positive $B$ is given by
eq.(43) with the last term replaced by $Bq^{2}$.

For the massive mode with the excitation energy $E_{+}$, the ${\bf
r}_{0}$-dependence is not important so that all of the expressions defined in
eq.(42) may be of unity. Then, the Hamiltonian expressing the Gaussian
amplitude fluctuation around the vortex lattice solution may be approximately
expressed as

$\displaystyle\delta{\cal H}_{\rm amp}$ $\displaystyle=$ $\displaystyle
N(0)\sum_{q,{\bf
k}_{\perp}}\biggl{(}\frac{{\tilde{V}}_{4}^{2}}{{\tilde{V}}_{6}}\biggl{(}1-{\overline{a}}^{\prime}_{0}-s_{4}\sqrt{1-{\overline{a}}^{\prime}_{0}}\biggr{)}$
(44) $\displaystyle+$ $\displaystyle
2|B|q^{2}+Cq^{4}\biggr{)}|\delta\Delta_{q}|^{2}$

when $B\leq 0$. This expression will be used in determining the melting line
through comparison between the free energy estimated from higher fields and
the corresponding one from lower fields.

As shown previously Moore ; RI591740 , the massless mode with the excitation
energy $E_{-}$ reduces to the purely phase fluctuation in $|{\bf
k}_{\perp}|\to 0$ limit and corresponds to the shear elastic mode of the
vortex lattice, where ${\bf k}_{\perp}=({\bf
r}_{0}\times{\hat{z}})/r_{H}^{2}$. By examining the ${\bf r}_{0}$-dependence
of the quantities defined in eq.(42) in details, the Hamiltonian expressing
the massless mode is found to become eq.(8).

## References

  * (1) A. A. Abrikosov, Sov. Phys. JETP 5, 1174 (1957).
  * (2) P. Fulde and R. A. Ferrell, Phys.Rev.135, A 550 (1964).
  * (3) A. I. Larkin and Yu. Ovchinnikov, Sov. Phys. JETP 20, 762 (1965).
  * (4) R. Ikeda, T. Ohmi, and T. Tsuneto, J. Phys. Soc. Jpn. 58, 1377 (1989).
  * (5) D. S. Fisher, M. P. A. Fisher, and D. A. Huse, Phys. Rev. B 43, 130 (1991).
  * (6) G. Blatter and V. B. Geshkenbein, in Superconductivity Vol.1 (ed. by K. H. Bennemann and J. B. Ketterson, Springer, 2008).
  * (7) A. Bianchi, R. Movshovich, C. Capan, P. G. Pagliuso, and J. L. Sarrao, Phys. Rev. Lett. 91, 187004 (2003).
  * (8) M. Kenzelmann, Rep. Prog. Phys. 80, 034501 (2017).
  * (9) S. Kasahara, Y. Sato, S. Licciardello, M. Culo, S. Arsenijevic, T. Ottenbros, T. Tominaga, J. Boker, I. Eremin, T. Shibauchi, J. Wosnitza, N. E. Hussey, and Y. Matsuda, Phys. Rev. Lett. 124, 107001 (2020).
  * (10) S. Kasahara, H. Suzuki, T. Machida, Y. Sato, Y. Ukai, H. Murayama, S. Suetsugu, Y. Kasahara, T. Shibauchi, T. Hanaguri, and Y. Matsuda, Phys. Rev. Lett. 127, 257001 (2021).
  * (11) F. Hardy, L. Doussoulin, T. Klein, M. He, A. Demuer, R. Willa, K. Willa, A.-A. Haghighirad, T. Wolf, M. Merz, C. Meingast, and C. Marcenat, Phys. Rev. Research 2, 033319 (2020).
  * (12) D. Y. Kim, S-Z. Lin, F. Weickert, M. Kenzelmann, E. D. Bauer, F. Ronning, J. D. Thompson, and R. Movshovich, Phys. Rev. X 6, 041059 (2016).
  * (13) H. Adachi and R. Ikeda, Phys. Rev. B 68, 184510 (2003).
  * (14) M. A. Moore, Phys. Rev. B 39, 136 (1989).
  * (15) S. Hikami, A. Fujita, and A. I. Larkin, Phys. Rev. B 44, 10400 (1991).
  * (16) R. Ikeda, Phys. Rev. B 76, 054517 (2007).
  * (17) G. Eilenberger, Phys. Rev. 164, 628 (1967).
  * (18) R.Ikeda, T. Ohmi, and T. Tsuneto, J. Phys. Soc. Jpn. 59, 1740 (1990).
  * (19) G. J. Ruggeri and D. J. Thouless, J. Phys. F6, 2063 (1976).
  * (20) K. Adachi and R. Ikeda, Phys. Rev. B 96, 184507 (2017).
  * (21) R.Ikeda, T. Ohmi, and T. Tsuneto, J. Phys. Soc. Jpn. 59, 1397 (1990).
  * (22) K. Maki, Phys. Rev. 148, 362 (1966).
  * (23) A. E. Koshelev, K.Willa, R.Willa, M.Smylie, J.-K.Bao, D.Y. Chung, M.G.Kanatzidis, W.-K. Kwok, and U. Welp, Phys. Rev. B 100, 094518 (2019).
  * (24) A. Houghton, R.A.Pelcovits, and A. Sudbo, Phys. Rev. B 40, 6763 (1989).
  * (25) K. Adachi and R. Ikeda, Phys. Rev. B 99, 064508 (2019).
  * (26) U. Klein, D. Rainer, and H. Shimahara, J. Low Temp. Phys. 118, 91 (2000).
  * (27) K. Yang and A. H. MacDonald, Phys. Rev. B 70, 094512 (2004).
  * (28) Y. Matsunaga, N. Hiasa, and R. Ikeda, Phys. Rev. B 78, 220508(R) (2008).
  * (29) S. Kasahara, T. Watashige, T. Hanaguri, Y.Kohsaka, T.Yamashita, Y.Shimoyama, Y.Mizukami, R.Endo, H.Ikeda, K.Aoyama, T.Terashima, S.Uji, T.Wolf, H.v.Lohneysenn, T.Shibauchi, and Y. Matsuda, Proc. Natl. Acad. Sci. USA 111, 16309 (2014).
  * (30) K. Adachi and R. Ikeda, J. Phys. Soc. Jpn. 84, 064712 (2015).
  * (31) D. F. Agterberg, Z. Zheng, and S. Mukherjee, Phys. Rev. Lett. 100, 017001 (2008).
  * (32) D. Nakashima and R. Ikeda, unpublished.
  * (33) N. Nunchot, D. Nakashima, and R. Ikeda, arXiv: 2203.11583.