Elliptical vortex and oblique vortex lattice in the FeSe superconductor based on the nematicity and mixed superconducting orders

The electronic nematic phase is characterized as an ordered state of matter with rotational symmetry breaking, and has been well studied in the quantum Hall system and the high-$T_c$ superconductors, regardless of cuprate or pnictide family. The nematic state in high-$T_c$ systems often relates to the structural transition or electronic instability in the normal phase. Nevertheless, the electronic states below the superconducting transition temperature is still an open question. With high-resolution scanning tunneling microscope measurements, direct observation of vortex core in FeSe thin films revealed the nematic superconducting state by Song \emph{et al}. Here, motivated by the experiment, we construct the extended Ginzburg-Landau free energy to describe the elliptical vortex, where a mixed \emph{s}-wave and \emph{d}-wave superconducting order is coupled to the nematic order. The nematic order induces the mixture of two superconducting orders and enhances the anisotropic interaction between the two superconducting orders, resulting in a symmetry breaking from $C_4$ to $C_2$. Consequently, the vortex cores are stretched into an elliptical shape. In the equilibrium state, the elliptical vortices assemble a lozenge-like vortex lattice, being well consistent with experimental results.

The electronic nematic phase is characterized as an ordered state of matter with rotational symmetry breaking, and has been well studied in the quantum Hall system and the high-Tc superconductors, regardless of cuprate or pnictide family. The nematic state in high-Tc systems often relates to the structural transition or electronic instability in the normal phase. Nevertheless, the electronic states below the superconducting transition temperature is still an open question. With high-resolution scanning tunneling microscope measurements, direct observation of vortex core in FeSe thin films revealed the nematic superconducting state by Song et al. Here, motivated by the experiment, we construct the extended Ginzburg-Landau free energy to describe the elliptical vortex, where a mixed s-wave and d -wave superconducting order is coupled to the nematic order. The nematic order induces the mixture of two superconducting orders and enhances the anisotropic interaction between the two superconducting orders, resulting in a symmetry breaking from C4 to C2. Consequently, the vortex cores are stretched into an elliptical shape. In the equilibrium state, the elliptical vortices assemble a lozenge-like vortex lattice, being well consistent with experimental results.

INTRODUCTION
In the newly discovered high superconducting transition temperature (T c ) iron-based family, the FeSe superconductors possess the simplest crystalline structure but attract much attention owing to multifarious physical properties [1][2][3][4] . The T c of bulk FeSe crystal is as low as 8 K, while it can be considerably enhanced to above 37 K under high-pressure 5 , electric field gating 6 , or insetting the intercalation layer 7 . Particularly, the monolayer FeSe on SrTiO 3 was observed a dramatically high T c above the liquid point of nitrogen 8 , which offers the possibility of breaking the record as those of cuprate family. The origin for the enhancement of T c is still an open question, while a common consensus has been proposed as an accompaniment to the modification of the Fermi surface. Therefore, studying on the electronic state of the FeSe system provides a perfect arena to understand the high-T c mechanism.
Different from the conventional superconductors, competing electronic orders such as unidirectional charge density wave and nematic order exist in both cuprate and iron-based superconductors. Among these, the nematic electronic order demonstrates a spontaneous symmetry breaking from C 4 to C 2 symmetry (the order parameters remain invariant under the inversion, the D 4h group can be viewed as C 4 ), which has been generally considered as a strong correlation with the fundamental unsolved electronic issue in Fe-based superconductors, especially in recent work on the FeSe system 9,10 . For the FeSe bulk crystals, the structural transition from tetragonal to orthorhombic occurs at T s = 90 K, while the anisotropy of the electronic structure is not a consequence of the lattice distortion, but a result of the microscopic mechanism such as spin fluctuation or orbital ordering. Researches on the nematic order in iron-based superconductors have generally supported the spin-fluctuation origin. However, because of the absence of long-range magnetic order in the FeSe system, orbital ordering is probably the origin for the electronic transition. Moreover, recent angle-resolved photoemission spectroscopy (ARPES) results showed the emergence of the nonequivalent energy shifts of xz/yz orbital bands below T s 10,11 , implying the orbital origin of the structural transition. Furthermore, similar to the nematic order, the superconducting pairing symmetry strongly relates to detailed electron-electron interaction. To be specific, the orbital order with inter-orbital electron-electron interactions would favor a sign-preserving s-wave pairing, while spin fluctuation with intra-orbital interaction for a sign-changing s ±wave or d-wave. Recent ARPES and Scanning Tunneling Microscopy (STM) results suggested the sign-changing pairing symmetry such as s ± -wave or d-wave in FeSe, implying that the magnetic fluctuations may still assist the superconducting pairing 12,13 . With the high-resolution STM measurement, the elliptical vortices have been directly observed at superconducting state on the FeSe bulk samples 14 , for which the extremely weak structure distortion (∼ 0.5%) can hardly induce such pronounced anisotropy, while the nematic order and superconducting order parameters are expected to play the important roles. Ginzburg-Landau (GL) theory firstly can offer a phenomenological way to investigate the vortices in superconductors with the s-wave symmetry. The GL theory itself can be derived exactly from the microscopic BCS theory 15 . By means of the Gorkov's derivation and symmetry analysis, the GL theory has been generalized into several pairing symmetries such as s + id [16][17][18][19][20] , p-wave 21 , and so on. Among these symmetry models the s+id, where the extended s ± -wave competes with the d-wave pairing order, is generally used to investigate the iron-based superconductors [22][23][24] .
In this work, we construct the GL type free energy which contains the nematic order, s-wave and d-wave superconducting orders with up to 4 th order interactions. The time-dependent GL (TDGL) equation is derived from the free energy to describe the FeSe system. By implementing the open boundary condition, our simulation reveals the configuration and dynamics of the elliptical vortex and the nematic order. With the periodical boundary condition, the oblique vortex lattice rather than a triangular one is found. Our simulation results have a good agreement with the previous experiment in the configuration of the single vortex and the vortex lattice 14 . The presence of the nematic order can break the symmetry from C 4 to C 2 and enhance the superconductivity. The symmetry allowed trilinear term will enhance the anisotropy of the superconducting order and induce the nearly degeneracy of s-wave and d-wave 25,26 .

RESULTS
High-resolution STM and scanning tunneling spectroscopy (STS) experiments provide the possibility for further investigation on the single vortex, for which the vortices configuration can be reconstructed as well studied in various superconductors [27][28][29][30][31][32] . In the previous work by Song et al. 14 , the vortices and vortex lattice in the FeSe superconductors were directly observed, and the vortex core was found in an elliptical shape, where the stretched direction is along one of the Fe-Fe bonds.
With the open boundary condition, the interplay between the anisotropic vortices and finite geometric region are investigated in the present work. Although one can hardly observe the evolution and dynamics of the vortices and nematic order in realistic experiments, the real-time simulation results can provide an approach to understand the motion of the vortices. By solving the TDGL equations, the results show that nematic order breaks the symmetry from C 4 to C 2 during the evolution. By using the periodic boundary condition, the vortex lattice is also investigated. Based on the simulation, the vortices favor an oblique lattice rather than the triangular lattice; this is due to the trade-off between the twofold symmetry of the repulsive interaction and the closet packing. The simulation results are consistent with the experimental data.

FINITE REGION AND VORTEX CONFIGURATION
Previous works suggest that the pairing symmetry is probably s ± wave or d wave, they both can be described by the addition of the isotropic and the anisotropic superconducting order. The isotropic order parameter is coupled to the anisotropic order parameter by the interaction 16,18,33 , where the two complex fields ψ s , ψ d stand for the s-wave component and d -wave component in the mixed superconducting order, and γ is the coupling constant. Π = (−i ∇ − e * A) is the gauge invariant derivative, ∇ is the del operator, A is the magnetic vector potential, and e * =2e is the charge of the superconducting charge-carriers, where e is the electron charge. This term is invariant under rotation of π/2, when taking the integration by parts, and thus, The above type of anisotropic interaction is used in our model. Besides the mixed superconducting order, a real field φ stands for the nematicity order, which competes with the mixed superconducting order. Because the higher order terms are negligible, the free energy which is up to 4 th order can be described as where α i , β i are the parameters describing the Landau phase transition and i = s, d and φ. Considering T s < T d < T φ , the α i submits to α s < α d < α φ . γ j (j = 1, 2, 3) is the coupling constant between s-wave and d-wave components, and λ k (k= 1-3) is the coupling constant between the superconducting order parameters and the nematic order. m * i (i = s, d) is the mass of the superconducting charge-carriers, and the microscopic electron pairing theory of superconductivity implies that m * i =2 m i , where m i is the electron mass. m φ represents the effective mass of the nematic order.
Instead of directly making the isotropic order parameter coupled to the nematic order, our model (Eq. 4) suggests that the nematic order triggers off the mixture of s-wave and d-wave components, and the anisotropic interaction between s-wave and d -wave components causes larger anisotropy. Meanwhile, the trilinear term λ 1 φ (ψ * s ψ d + c.c.) will enhance the anisotropy. Besides, according to the previous mean-field analysis, this term may induce nearly degeneracy of s-wave and d-wave, which was also supported by the spin-fluctuation model 25,26 . This is different from the p-type 122-system, where the degeneracy of s-wave and d-wave is due to the close critical temperatures 34 .
In the following simulation, the phase transition parameters in the Eq. 4 are set to be α s = 1.0, α d = 1.5, α φ = 2.0, β i = 1, i = s, d, and φ, which are based on the superconductivity and nematicity transition temperature 9,12,35 . The coupling constant of the anisotropic interaction is set as γ 3 = 0.2, while other coupling constants are set to be 0.4 based on the consideration of convergency. λ 2 , λ 3 are negative. The effective masses are set to be m s = 1, m d = 2, m φ = 4.
The initial condition is quite important for convergency of the non-linear partial differential equations, though different initial conditions will arrive at the same stable state in this type partial differential equations. The initial states of the complex order parameters are set to be proportional to For the s-wave component, the initial state is ψ 0 , while it is relatively small for the d-wave component. The initial state of the nematic order is set to be 0.5, and all the vector potentials are set to be 0.

EVOLUTION OF THE VORTICES AND THE NEMATIC ORDER
The vortex core can be visualized by atomically resolved STM measurements. Moreover, by applying the specific periodical magnetic field, it is possible to view the motion of the vortex 36 . However, it is still difficult to conduct real-time observation on the evolution and dynamics of the vortices. Time-dependent simulation can be applied to simulate how the vortices generate from the boundary and how they move and interact with each other. Through the real-time simulation, deep understanding of the vortex dynamics can be achieved and further applications can be simulated, such as modification of the sample to enhance the critical current.
Different from the other type-II superconductors, the elliptical vortices in FeSe sample is related to the C 4 → C 2 symmetry breaking, where the nematic order plays an important role on enhancing the symmetry breaking in the superconducting state 37 . With the real-time simulation, the transition from C 4 symmetry to C 2 symmetry is revealed (see Supplementary videos).
The real-time evolution of the elliptical vortices is shown in Fig. 1. At the beginning, the pattern of the s-wave component shown in Fig. 1a.1 is quite similar to that in typical type-II superconductors [38][39][40][41] , the magnetic field penetrates into the sample from the edge. Due to the Bean-Livingston barrier 42,43 , the vortices cannot immediately get into the sample. The interaction of intra-vortices is qualitatively repulsive and the vortices could only locate along the edges.
The enhanced magnetic field will force the vortices to penetrate into the sample. As a result, Fig. 1a.2 demonstrates the arrangement of four vortices. The system will gradually achieve its minimum free energy by rearranging the vortices. However, with the nematic order competing with the mixed superconducting order, the situation is quite different. Fig. 1a.3-6 depict the intermediate stage, where the C 4 → C 2 symmetry breaking happens, because the nematic order mixes the two superconducting order and enhances the anisotropic interaction. Meanwhile, the interaction terms λ 1 and γ 1,2 compete with each other, λ 1 favors a large separation between s-wave and d-wave, while the other one favors small separation. Fig. 1a.7 demonstrates the equilibrium state of this system, the two vortices are both elliptical and repulse with each other in short range.
Different from the s-wave component, the initial state of the d-wave component is relatively small and then induced to be anisotropic by the nematic order. The vortices of d-wave component are slightly less eccentric due to its fourfold tendency. Apart from the superconducting order, the nematic order shown in Fig. 1c.1-7 exhibits strong C 2 symmetry at first. It then has fourfold symmetry due to the interaction with s-wave component. The final state shown in Fig.  1c.7 has C 2 symmetry and is less eccentric than the superconducting order.

VORTEX LATTICE
Vortices arrange into the vortex lattice in an infinite region. Though it is impossible to simulate the complicated TDGL equations in an infinite region, simulation on the single unit cell with periodic boundary condition can make the investigation of the vortex lattice possible [44][45][46] . By extending the unit cell according to the periodicity, the vortex lattice can be recovered.
Defining the ratio of the side lengths of the rectangular unit cell as r, r = √ 3 corresponds to triangular lattice which is typical for most type-II superconductors, as shown in Fig. 2a. Early study reported that the oblique vortex lattice 16,44,45,47 , r < √ 3, for s + id model costs less energy than the triangular one which is ascribed to the fourfold symmetry of the system. To be specific, the system tends to preserve the fourfold symmetry while the closet packing between the vortices leads the vortex lattice to be triangular. As a result, the vortex lattice favors the oblique one by making the trade-off between preserving the fourfold symmetry and the triangular lattice. While in our case, the nematic order breaks the C 4 symmetry to C 2 , and thus the system finds the balance between the twofold symmetry and the triangular lattice. Therefore, it favors the lattice with r > √ 3.
In the simulation, the normalized magnetic field is set as 4π, which allows two vortices in the one unit cell. The area of the unit cell is set to be 16λ 2 . By varying the ratio r, the minimum energy density f = F/L x /L y is achieved at r = 2.81, as shown in Fig. 2. It is hard to realize the real-time detection in experiments and compare the vortex dynamics with the simulation, but the equilibrium state where the vortices form a stable vortex lattice can be compared with the simulation results. Based on the parameters provided above, the simulated oblique vortex lattice r = 2.81 is in agreement with the previous result r ∼ 2.80 in the experiment where many vortices are observed in FeSe under applied magnetic field of 8 T 14 .

ANALYTICAL TREATMENT OF THE ANISOTROPY
According to the previous simulation results, the anisotropy of the interaction serves to the elliptical shape of the vortices. Since the London penetration depth is considerably large than the coherence length λ ξ, the coupling to the electromagnetic field can be neglected. Given that the variation of the nematic order is small, and the vortices experience a uniform nematic order away from the origin, thus nematic order is set to be a stationary field, φ = φ 0 , and satisfy φ 0 → −φ 0 under the rotation of π/2. The equations of the superconducting order are, Where γ s,d = 2 2m s,d . A direct observation on Eq. 5 reveals that λ 2 , λ 3 and φ 2 0 will change the critical point of the phase transition in this system, as the minimum point of the potential is at ψ s,min = ± α s − λ2 The coefficients of the term ψ s and ψ d are non-zero and different for most situations. Because λ 2 , λ 3 are negative, the presence of the nematic order will enhance the superconductivity 48 .
The interesting question is that the λ 1 term will break the rotation symmetry. When rotating π/2, the λ 1 term picks up a different sign compared with other terms. Such term enhances the anisotropy of the vortices.
Polynomial terms in Eq. 5 do not contribute to the anisotropy but the gradient terms and the nematic order will. Without losing universality, we set γ s = γ d = γ. Eq. 5 can be reformulated as, where ψ = ψ s +ψ d , P is the polynomial of ψ s and ψ d . Ignoring the nematic order φ 0 , the solution is actually a elliptical vortex, which can be seen by transforming the gradient terms to a Laplacian under the coordinate transformation, Due to Eq. 7, the coordinate is elongated along y direction and contracted along x direction as shown in Fig. 3. Therefore, the symmetry is broken into C 2 once returning to the original coordinate. However, the nematic order obeys φ 0 → −φ 0 under rotation of π/2, it is impossible to view x-direction and ydirection equivalently, because the nematic order offers a angle dependent term. To capture the feature, φ 0 is set to be k(x 2 − y 2 ), where k is a constant. Fig. 4 shows that turning on the nematic order makes the elliptical vortex more anisotropic.

COMPARISON WITH OTHER SYSTEMS
The anisotropic electronic structure is normal in some high-T c superconductors. For instance, in the cuprate family, YBa 2 Cu 3 O 7−δ (YBCO) has a tetragonal to orthorhombic phase transition at the under-doped level of oxygen, resulting in a symmetry breaking from D 4h to D 2h along the Cu-O chains 49 . Here, the structure transition temperature is considerably higher than the superconducting transition temperature T c . The previous works on the YBCO have constructed the GL free energy which obeys D 2h symmetry and found the elliptical vortex along the b-axis 50,51 . In the present FeSe system, orthorhombic phase happens after a structural phase transition at 90 K. However, it is argued that the small crystalline distortion (∼ 0.5%) itself cannot lead to such large anisotropy in electronic structure 14,52 . Therefore, the nematic order may originate from the strongly interacting fermion system, where the small symmetry breaking term will be promoted to finite magnitude when lowering the energy scale, say, temperature. Thus, different from the explicit symmetry breaking treatment in YBCO, the symmetry breaking nematic order in FeSe is temperature dependent and should be treated as a competing order which competes with the D 4h superconducting order 24 . On the other hand, the elliptical vortex core is along the Fe-Fe bonds, being considerably different from that along the b−axis in YBCO.
In some iron-based superconductors, the superconducting nematic transition temperature was also found below the superconducting transition temperature, such as p-type 122-system iron pnictide superconductors 34 . The nematic superconducting state in this superconductor is induced by the small symmetry breaking term in the normal phase which qualitatively differs from the ordinary nematicity observed in the orthorhombic structural phase. The superconducting state nematicity occurs just after the onset of superconducting transition and reveals an anisotropy shifted by π/4 from those of the normal nematic state 34 . Meanwhile, the superconducting transition temperature resembles a low energy scale, hereby, the lower nematic transition temperature may correspond to an infrared phenomenon of the strongly interacting fermion system. However, one can hardly identify which order occurs first. When the superconducting transition temperature T c is lower than the nematic transition temperature T n , it can be expected that the nematic order still exists in superconducting state and strongly interacts with the superconducting orders, which are shown in Fig. 5a. For the case of T c > T n , there is an intermediate region T c > T > T n , where the nematic order is absent in the superconducting state, and the profile of the vortex core is shown in Fig. 5b. Without the nematic order, the mixture of s-wave and d-wave makes the vortex slightly anisotropic. However, with the nematic order occurring below T n , the anisotropy of the superconducting vortex is enhanced. In the simulation, the high order correction term is added into the free energy.

METHODS
We begin by constructing the GL type free energy. Intuitively, the nematic order competes with the mixed superconducting order parameters, where isotropic s-wave and anisotropic d -wave components are considered. We mainly focused on the 2-dimensional (2D) geometry due to the quasi-2D feature for the Fe-based superconductors, and the 2D GL free energy can mostly capture the ingredients in the high-T c superconductors. By taking the variation of the free energy, the TDGL equations are derived. The TDGL was constructed to investigate the dynamics of the vortices in the dirty limit where the penetration depth λ is greatly larger than the coherence length ξ 53 . For the FeSe case, due to λ ∼ 500 nm and ξ ∼ 5 nm 54,56? , the dirty limit construction of TDGL is valid. Numerical results of the TDGL equations with the open boundary condition, given by the finite element method, can reveal the shape, configuration and dynamic properties of the vortices. The numerical solutions with periodic boundary condition provide the vortex lattice in the equilibrium state. Our theoretical calculation results show the elliptical vortices and oblique vortex lattice, which are in agreement with the experiment 14 . In the following formalism, the a− and b−axis, or x− and y−axis are defined along either of the Fe-Fe bond directions as shown in Fig. 6, for which the directions of the Fe-Fe bonds keep the symmetry of the nematicity.

FREE ENERGY
Competing order such as nematicity can strongly interact with the superconducting order parameters. Chowdhury et al. investigated the anisotropic interplay between the competing order and the single isotropic superconducting order 23 , and they argued that the different effective masses which are induced by the anisotropic interaction could lead to the anisotropic vortex. However, the pairing symmetry in FeSe is sign-changing s ± -wave or d-wave suggested by the recent experiments 12,13 . Meanwhile, the s + id model could suitably describe the iron-based superconductors 22 . Here, the mixed superconducting order together with the anisotropic interaction can also cause an anisotropic vortex core, where the isotropic s-wave order parameter interacts with the anisotropic d -wave order parameter. The existence of nematic order will mixed the two superconducting order and significantly enhance the anisotropy. Thus, the nematic order can enhance the small anisotropic interaction between the superconducting orders to form an extremely elliptical vortex.
Previous work investigated the special trilinear term thoroughly. It turns out if the coupling constant λ 1 = 0 or very small, the system may favor s + id symmetry, but for large nematic fluctuation, the intermediate state has s + d symmetry 26 character. The free energy (Eq. 4), including both self-energy and interaction energy, remains invariant under the mirror reflection and the rotation of π/2.
However, the γ 3 term in the interaction causes different effective masses along the two directions, and consequently, results in the anisotropic vortex cores and affects the arrangement of the vortices, namely, the vortex lattice. Based on the consideration of symmetry, up to 4 th order, it is impossible to turn on the direct interaction between the nematic order and the anisotropic gradient term. Nevertheless, the nematic order can tune the stationary part of the superconducting orders in the free energy and let them mixed, thus the anisotropic interaction γ 3 enhanced. The λ 1 term will also enhance the anisotropy as explained in Section Analytical Treatment of the Anisotropy. Because the λ 2 , λ 3 are negative, the presence of the nematic order will enhance the superconductivity, such as increasing the transition temperature.

TIME-DEPENDENT GINZBURG-LANDAU EQUATION
The TDGL equation can be obtained by taking the variation of the free energy as follow, where D i (i = s, d, φ) is the phenomenological diffusion coefficients, and Φ is the scalar potential of the electromagnetic field.
In the open boundary conditions, the superconductivity-vacuum boundary can be obtained directly from the variation of the free energy as wherex,ŷ are the unit vector along x and y directions, respectively, and l x,i = 1+ mi m sd and l y,i = 1− mi m sd correspond to the different effective mass along x and y directions, respectively. Thus, the dynamics property of the superconducting order is different along the x and y axis.
To solve the TDGL equations numerically, the complicated TDGL equations are normalized by introducing, where the new quantities are labeled by prime, the spatial and temporal coordinates are scaled according to the λ and the ξ = √ 2msαs , the GL parameter is defined as κ = λ/ξ and κ 1 for the FeSe system 54,56? . The dimensionless form contains the gauge invariant derivative, Π = − i κ ∇ − A. The TDGL is invariant under the gauge transformation. Given an arbitrary function χ (x, y, z, t), and introduce the gauge transformation as,ψ Because of the extra degree of freedom, the gauge should be fixed to achieve the definite equations. For the sake of the simplicity, the scalar potential Φ can be eliminated by choosing the London gauge, let, thus the TDGL equations are no longer dependent on the scalar potential Φ, and the electric field is E = − ∂A ∂t . The former procedure derives the TDGL equations and corresponding open boundary conditions based on the free energy. To implement the finite element method, the complex order parameters are decomposed into the real and imaginary part; the vector potentials are decomposed into x, y components. For the consequences of the non-linear feature, the mesh of the region is adaptively refined to achieve high accuracy. Solving TDGL is minimizing the total energy of the system, and the stable state will be reached after hundreds to thousands of the normalized time.

PERIODIC BOUNDARY CONDITION
The open boundary condition describes the superconductor-vacuum boundary straightforward, therefore, it is useful to investigate the finite size solution. When coming to the infinite size, the periodic boundary condition should be introduced on each unit cell. The complex order parameter and the vector potential should be modified from one unit cell to another which can ensure the gauge invariance 57 . Two lattice vectors are used to characterize the lattice, namely, t 1 and t 1 . The complex order parameters will pick up a phase while additional term should be added to the vector potential from one unit cell to another, namely, where, g k = − 1 2 (t k × Bk 3 ) · x, k = 1 and 2, and B is quantized by B = 2πn κ|Ω| , in which |Ω| is the area of the unit cell and n is integer.
The explicit form for a rectangular unit cell is, ψ i (L x , y) = ψ i (0, y) e iφy 2Ly A x (L x , y) = A x (0, y) A y (L x , y) = A y (0, y) + φ where, i = s and d, φ = 2nπ is the reduced vortex flux, L x and L y characterize the size of the unit cell. The variation of the vector potential is neglected, due to κ 1.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.