Intertwined Orders and Electronic Structure in Superconducting Vortex Halos

The high-temperature superconductivity (SC) of cuprates is marked by the many ordering tendencies that either compete or coexist with the superconducting phase itself Keimer et al. (2015); Fradkin et al. (2015); Martin et al. (1990); Daou et al. (2009); Varma (2020). For example, charge density wave (CDW), either unidirectional or bidirectional, have been observed in all families of hole-doped high-temperature superconductors (HTSCs) Tranquada et al. (1995); Kivelson et al. (2003); Comin and Damascelli (2016). The extensively studied stripe order in cuprates correspond to a unidirectional CDW, often accompanied by a spin density wave (SDW). The CDW order, which is intimately related to the intriguing pseudogap phase of cuprates, is short-ranged at zero field Blackburn et al. (2013); Ghiringhelli et al. (2012); Blanco-Canosa et al. (2014); Tabis et al. (2017). Both the strength and correlation length of the CDW is enhanced by the magnetic field Wu et al. (2011, 2013); Zhou et al. (2017). On the other hand, the onset of SC causes the reduction of the CDW amplitude, suggesting that these two orders are strongly intertwined with each other Ghiringhelli et al. (2012); Chang et al. (2012). The many unusual features of charge-density modulations in the pseudogap phase, such as the doping dependence of the modulation period, suggest that the CDW is a subsidiary order which results from a more fundamental pair-density wave (PDW) ordering.

A PDW is a superconducting state in which the pairing order parameter varies periodically in space Agterberg et al. (2020). The idea of a modulated pairing density was first introduced by Fulde and Ferrell and by Larkin and Ovchinnikov (FFLO) in a BCS model to overcome the Pauli limiting effect of a magnetic field Fulde and Ferrell (1964); Larkin (1965). Recently, PDW states without any explicit breaking of the time-reversal symmetry have been proposed to exist in cuprates, especially in connection with the pseudogap phases Chen et al. (2004); Berg et al. (2009a); Lee (2014). In particular, the PDW order and its partial melting could lead to a variety of vestigial states including the CDW, possibly coexisting with an SDW mentioned above, and an unusual charge-$4e$ superconductivity Berg et al. (2009b), among others. The rich PDW phenomenology seems to provide a natural explanation for the complexity of cuprate HTSC. Moreover, the scenario of decoupled layers of orthogonal planar PDW orders could account for the observed huge anisotropy of resistivity well above the nominal SC transition temperatures in La_2-xBa_xCuO₄ at the 1/8 hole doping Li et al. (2007); Tranquada et al. (2008); Berg et al. (2007).

Direct imaging of PDW orders has also been reported in hole doped Bi₂Sr₂CaCu₂O_8+x (Bi2212) using atomic-resolution superconducting STM tips Du et al. (2020). The amplitude of these modulations is characterized by a eight-unit-cell periodicity or wavevectors $\mathbf{Q}_{x}=\left(1/8,0\right)(2\pi/a)$ and $\mathbf{Q}_{y}=\left(0,1/8\right)(2\pi/a)$, where $a$ is the lattice constant. Simultaneous measurement of local density of states found electronic modulations with wavevectors $\mathbf{Q}_{x,y}$ and $2\mathbf{Q}_{x,y}$ Du et al. (2020), which are consistent with the picture of a PDW coexisting with a uniform $d$-wave superconductivity Agterberg and Garaud (2015a). The nature of this SC state at zero magnetic field and its relationship with the PDW order, however, remains to be resolved. Given the complex nature of the cuprate superconductors with several energetically nearly degenerate ordered states, extrinsic effects such as impurity and disorder likely play an important role in the stabilization of this intriguing state with mixed SC and PDW ordering.

As the pairing order parameter is suppressed near a vortex core, magnetic-field induced vortices offer a fruitful platform to further investigate the subtle interplay between PDW and superconducting condensate. Indeed, for quite some time, a number of STM experiments have reported a nonuniform charge density and spectra inside the halo region surrounding the vortex core Hoffman et al. (2002); Matsuba et al. (2003, 2007); Hanaguri et al. (2009); Hamidian et al. (2015); Machida et al. (2016). A checkerboard pattern with a $4a\times 4a$ unit cell were observed in almost all samples. This checkerboard charge density modulation can also be viewed as a bi-directional CDW with wavevectors $2\mathbf{Q}_{x}$ and $2\mathbf{Q}_{y}$, which are twice that of a fundamental PDW discussed above Du et al. (2020); Edkins et al. (2019). What is also striking is the bipartite electronic structure at the vortex core Machida et al. (2016). At low energy quasiparticle interference (QPI) dominates in the conductance while at higher energy, the energy-independent CDW patterns emerge.

Recently this vortex halo has been further studied and proposed to be a bidirectional PDW by Edkins et al. Edkins et al. (2019). The discrete Fourier transform of the conductance map of the halo reveals a modulation of the charge density $N(\mathbf{r})$ with peaks at wavevector $\mathbf{Q}_{x,y}$, in addition to the previously observed $2\mathbf{Q}_{x,y}$. Phenomenologically, this again suggests the coexistence of uniform SC condensate $\Delta_{\rm SC}$ and the PDW order $\Delta_{\mathbf{Q}}$ with wavevectors $\mathbf{Q}=\mathbf{Q}_{x}$ and $\mathbf{Q}_{y}$. In this scenario, a CDW order could result from couplings $\Delta^{\,}_{\mathbf{Q}}\Delta^{*}_{\rm SC}$ as well as $\Delta^{\,}_{\mathbf{Q}}\Delta^{*}_{-\mathbf{Q}}$, hence giving rise to charge modulations with a wavevector $\mathbf{Q}$ and $2\mathbf{Q}$, respectively Agterberg et al. (2020).

The recent STM experiments also reinvigorate a long-standing puzzle of cuprates. Due to the nodes of $d$-wave pairing order, it is expected and confirmed by theoretical calculations Wang and MacDonald (1995); Franz and Tešanović (1998) that a zero-biased conductance peak (ZBCP) not should be present at an isolated vortex center. Intriguingly, however, this peak has so far not been detected in most experiments. Instead, a sub-gap structure was reported Machida et al. (2016). Several scenarios, e.g. the presence of a concomitant antiferromagnetic or SDW order Zhu and Ting (2001); Franz et al. (2002) or an induced hidden $d_{xy}$ pairing  Franz and Tešanović (1998) in the halo region, have been proposed to explain this discrepancy. A consensus has yet to be reached among researchers.

The issue of the vortex-induced ZBCP is further complicated by a recent experiment showing that, in contrast to sub-gap structure under higher magnetic field at $B\sim 3$T, there is a ZBCP at the vortex core under extremely low magnetic field ($B\sim 0.16$ T) in an over-doped Bi2212 sample Gazdić et al. (2021). In addition, at a hole doping of around $\delta=0.2$, the conductance map and spectra in the halo region of a $3$T vortex are shown to be very similar to that due to a checkerboard charge modulation in under-doped samples reported by others. However, the wavevectors associated with the checkerboard seem to disperse within an energy range of 3 to 10 meV. This result seems to be related with the aforementioned bipartite structure found by Machida et al. Machida et al. (2016) for an optimally doped sample under 11.25 T magnetic field, yet with very different dopant density and field strength.

All these recent results posed several new interesting questions. Is the ZBCP only present in over-doped samples? Or is it suppressed by the presence of PDW or checkerboard state? Is the bipartite structure inside a vortex an intrinsic property of the PDW state? What role does the QPI play in the formation of checkerboard halo states? What is the mechanism for the inducement of the checkerboard pattern inside a vortex? To answer these questions, one must examine the vortex structure as a function of magnetic field and dopant density. In particular, under and optimal-doped samples may have very different structures from that of the over-doped ones.

The effects of PDW order on the vortex and the associated electronic structures have also been examined in several recent theoretical studies Agterberg and Garaud (2015b); Dai et al. (2018); Wang et al. (2018a). In these approaches, the order parameter fields of a vortex are obtained from phenomenological Ginzburg-Landau free energy functionals. A quadratic lattice fermion model with input from the order-parameter solutions is diagonalized to study the electronic properties, e.g. local density of states, of the vortex. Such semi-empirical approaches can provide qualitatively experimental signatures due to different ordering scenarios, such as whether the checkerboard pattens is driven by a primary PDW order, or the other way around. Yet, a fully self-consistent theoretical modeling of vortex halos, especially one based on well defined microscopic models, is still lacking. Moreover, electron correlation effects, which are known to be important for cuprate materials, cannot be properly included in the phenomenological approaches.

In this paper, we present a comprehensive study on the structure of vortex halos based on large-scale self-consistent calculations of the well-studied $t$-$J$ model with an additional frustrated second-neighbor $t^{\prime}$ hopping. It is worth noting that several advanced many-body techniques, ranging from unrestricted mean-field method and variational Monte Carlo (VMC) simulation to density matrix renormalization group (DMRG) and tensor-network methods, have been applied to study the intertwined orders in strongly correlated electron systems Loder et al. (2011); Himeda et al. (2002); Berg et al. (2010); Dodaro et al. (2017); Jiang et al. (2018); Zheng et al. (2017); Corboz et al. (2014). In particular, PDW order has been demonstrated in the Hubbard and $t$-$J$ models, as well as their variants. While sophisticated methods, such as DMRG or infinite projected entangled paired states (iPEPS), generally can give more accurate results, and numerically exact results for 1D models, so far they can only be applied to relatively small 2D systems, and hence are not feasible for studying large-scale inhomogeneous states such as superconducting vortices.

On the other hand, renormalized mean-field theory (RMFT) has proven an efficient method for solving complex structures in $t$-$J$ type correlated electron models Garg et al. (2008); Yang et al. (2009); Christensen et al. (2011); Tu and Lee (2016); Banerjee et al. (2018). As in most mean-field type approaches, the many-body Hamiltonian is reduced to an effective single-particle problem which is to be solved self-consistently. Yet, unlike standard Hartree-Fock type approximations, the strong electron correlation effects in such systems can be properly captured by the Gutzwiller approach employed in the RMFT Gutzwiller (1963); Himeda and Ogata (1999); Ogata and Himeda (2003). Importantly, fairly quantitative agreement with STM and ARPES experiments on cuprates have been obtained from RMFT calculations of the frustrated $t$-$t^{\prime}$-$J$ model Choubey et al. (2017a, 2020); Wang et al. (2021); Tu and Lee (2019); Tsuchiura et al. (2003). The efficiency of our RMFT calculations is further improved significantly by incorporating the kernel polynomial method (KPM) Wang et al. (2018b); Weiße et al. (2006) for solving the renormalized tight-binding Hamiltonians, which has to be repeated up to thousand of times in the self-consistency calculation. By applying the RMFT method to the the $t$-$t^{\prime}$-$J$ model on square-lattice of up to $10^{4}$ sites, we obtain detailed vortex solutions where up to $10^{5}$ variational parameters are solved self-consistently.

Our large-scale RMFT calculations find several vortex solutions with nearly degenerate energies. This result again underscores the complex nature of the superconducting phase, and is reminiscent of previous works showing multiple density-wave states whose energies are almost degenerate with that of uniform $d$-wave SC state Tu and Lee (2016); Choubey et al. (2017a); Zheng et al. (2017); Corboz et al. (2014). In addition to conventional vortices in a uniform $d$-wave background, of particular interest is a self-consistent vortex solution where a bidirectional PDW with a period $8a$ coexists with a charge-density modulation of the same period, which is twice the period of the charge-density wave order resulting from the second harmonic of a unidirectional PDW in the absence of a vortex. This solution, which is obtained without invoking special features of the Fermi surface, is consistent with recent STM experiment Edkins et al. (2019).

By carrying out systematic study of the field and doping dependences of these vortex solutions, a coherent and consistent picture is provided for recent STM experiments and the issues of ZBCP. We show that, due to the strong correlation and the orbital $d_{x^{2}-y^{2}}$ symmetry, the conductance spectra of STM have sub-gap structures instead of a ZBCP in the under- and optimal-doped regimes even for the plain vortex state. Moreover, the vortex solution with a bipartite checkerboard patterns is shown to result from coexisting bi-directional PDW, CDW, and SDW orders. Such vortices with intertwined orders exhibit energy independent wave vectors at higher energy but with QPI dispersion at low energy. It also shows a larger $s$+$s^{\prime}$ form factor than that of the $d$-wave symmetry. Finally, the mechanism for the emergence of the checkerboard state inside a vortex is discussed.

The hole-doped CuO₂ plane of copper oxides can be well described by the generalized $t$-$J$ model, which corresponds to the strong coupling limit of the Hubbard model. Its Hamiltonian reads

where $\hat{c}^{\dagger}_{i\sigma}$ ($\hat{c}^{\,}_{i\sigma}$) is the creation (annihilation) operator of an electron at site-$i$ with spin $\sigma$, and $\hat{\mathbf{S}}_{i}=\hat{c}^{\dagger}_{i\sigma}\mathbf{\sigma}_{\sigma,\sigma^{\prime}}\hat{c}^{\,}_{i\sigma^{\prime}}$ is the electron spin operator at site-$i$, $t_{ij}$ denotes the electron transfer integral between the Cu $d_{x^{2}-y^{2}}$ orbitals at sites $i$ and $j$. In the $t$-$t^{\prime}$-$J$ model to be studied in this work, only electron hopping $t$ between nearest-neighbor and $t^{\prime}$ between next-nearest-neighbor pairs are included in the first term. The notation $\langle ij\rangle$ in the second term denotes the nearest neighbors. The effect of the strong on-site Coulomb repulsion $U\to\infty$ is accounted for by the Gutzwiller projector $\hat{P}_{G}=\prod_{i}\left(1-\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}\right)$, which eliminates all doubly-occupied orbitals; here $\hat{n}_{i\sigma}=\hat{c}^{\dagger}_{i\sigma}\hat{c}^{\,}_{i\sigma}$ is the electron number operator at site-$i$. The residual charge fluctuations between nearest-neighboring pairs $\langle ij\rangle$ leads to the antiferromagnetic superexchange $J\sim t^{2}/U$ in the second term above. In the presence of a magnetic field $\mathbf{B}$, its effect is accounted for by the Peierls substitution $t_{ij}=te^{i\phi_{ij}}$, where $\phi_{ij}=\frac{-\pi}{\Phi_{0}}\int_{\mathbf{r}_{i}}^{\mathbf{r}_{j}}{\mathbf{A}}(\mathbf{r}_{i})\cdot d\mathbf{r}_{i}$, where $\Phi_{0}=\frac{hc}{2e}$ is the superconducting flux quantum, and $\mathbf{A}$ is the corresponding vector potential. In our implementation, a Landau gauge $\mathbf{A}=B(0,x)$ is used.

For large-scale calculations of inhomogeneous states such as those with vortices, the above $t$-$t^{\prime}$-$J$ model is solved using the RMFT method Garg et al. (2008); Yang et al. (2009); Christensen et al. (2011); Tu and Lee (2016); Banerjee et al. (2018); Choubey et al. (2017a); Tu and Lee (2019); Choubey et al. (2020); Wang et al. (2021); Tsuchiura et al. (2003). In this approach, the strong electron correlation, as encapsulated by the projector $\hat{P}_{G}$, is treated by the Gutzwiller approximation, giving rise to the following renormalized Hamiltonian

Here $g^{t}_{ij}$ and $g^{s}_{ij}$ are Gutzwiller factors that account for the renormalization of the hopping amplitude and the exchange interaction, respectively. The exchange interaction $J$, which leads to the SC pairing O’Mahony et al. (2022) and other intertwined orders, is treated by the Hartree-Fock-Bogoliubov mean-field decoupling. In this work, we consider the following mean-field order-parameters: (i) the on-site hole density

(ii) the local magnetization

(ii) the bond-order between neighboring Cu sites

and (iii) the local electron pairing field on nearest-neighbor bonds ($\bar{\sigma}=-\sigma$)

Here the superscript $v$ is a reminder that these parameters, while each related to physical observables, are not directly measurable quantities; a proper renormalization factors have to be included for experimental comparisons. The avereage $\langle\hat{O}\rangle=\langle\Psi_{0}|\hat{O}|\Psi_{0}\rangle$ denotes expectation value of operator $\hat{O}$ with respect to a variational Slater wave function $|\Psi_{0}\rangle$. The Slater state is constructed from eigenstates of an effective tight-binding Bogoliubov-de Gennes (TB-BdG) Hamiltonian

where $\mathcal{T}_{ij\sigma}$, $\mathcal{D}_{ij\sigma}$, and $\mu_{i\sigma}$ are effective hopping, pairing, and on-site energy parameters, respectively, defined as the derivatives of the variational energy of the renormalized Hamiltonian with respect to the local order parameters:

Here $W=\langle\Psi_{0}|\hat{H}_{\rm renorm}|\Psi_{0}\rangle$ supplemented by Lagrangian multipliers for enforcing the wave function normalization and the conservation of electron number, and $n_{i\sigma}=\langle\hat{n}_{i\sigma}\rangle$ is the local electron number; see Appendix A for details. These parameters in turn depend the local order-parameters listed in Eqs. (3)–(6), that are determined self-consistently. In order to allow for spatial inhomogeneity, these site and bond-dependent parameters, with a total number of the order of $10^{5}$, are treated as independent and optimized in real-space RMFT calculations; details of the method are presented in Appendix A.

In the RMFT calculations, the numerous local mean-field parameters are solved self-consistently through iterations. For systems with up to $N=10^{4}\sim 10^{5}$ sites, each iteration requires solving a large $2N\times 2N$ tight-binding matrix. For unrestricted optimizations, a number of $1000\sim 3000$ iterations are routinely required to reach satisfactory convergence. Although the TB-BdG Hamiltonian can be exactly solved by the exact diagonalization (ED), the poor $\mathcal{O}(N^{3})$ scaling of ED renders it infeasible for large-scale calculations. Here, instead, we incorporate the kernel polynomial method (KPM) Weiße et al. (2006); Barros and Kato (2013); Wang et al. (2018b), which provides a linear-scaling $\mathcal{O}(N)$ approach to electronic structure problems, into the RMFT framework.

In this approach, the total or local DOS of the system is expanded as a Chebyshev polynomial series whose coefficients are efficiently computed through matrix-vector products, which can be made linear-scaling for large sparse TB matrices. The various correlation functions $\langle c^{\dagger}_{i\sigma}c^{\,}_{j\sigma^{\prime}}\rangle$ and $\langle c^{\,}_{i\sigma}c^{\,}_{j\sigma^{\prime}}\rangle$ required for the computation of local order parameters are obtained through automatic differentiation. We note in passing that KPM has been used within the conventional Hartree-Fock-Bogoliubov mean-field to study large-scale structures of SC states Covaci et al. (2010); Nagai et al. (2012); Berthod (2016). Yet, the introduction of the Gutzwiller renormalization factors significantly increases the computational complexity of RMFT. In our implementation, the efficiency of KPM is further enhanced by employing the probing method Silver and Röder (1994); Tang and Saad for estimating the trace of large matrices and the utilization of GPU programming; more details can be found in Appendix B.

The KPM-RMFT method discussed above is mainly used to obtain the various local parameters in a vortex structure, a calculation which often requires up to thousands repeated solutions of the TB-BdG Hamiltonians in order to reach self-consistency. Once these local orders are obtained, a hybrid approach is employed to compute the corresponding electronic structures such as local density of states. We combine efficient exact diagonalization (ED) package Tomov et al. (2010) with the supercell method Zhu et al. (2002); Schmid et al. (2010) to obtain the exact eigenfunctions of the converged TB-BdG equation. We note that even though large-scale ED is time-consuming, this is a one-shot calculation using the converged order-parameters. The supercell method allows us to further reduce the finite-size effect and enhance the momentum resolution.

The electron Green’s function, which is directly related to several experimental measurements, can be computed from the exact eigenfunctions of the renormalized TB-BdG Hamiltonian. In particular, the local density of states (LDOS) at the $i$-th site corresponds to the diagonal element of the lattice Green’s function

For better comparisons with experiments, one can take into account the effects of electron orbitals using the continuum Green’s functions Choubey et al. (2014); Kreisel et al. (2015) defined as

Here $W_{i}(\mathbf{r})$ is the Wannier function centered at site $i$. The corresponding continuum LDOS, which can be directly measured by the tunneling conductance, is

In practical calculations, a broadening factor $\Gamma$ is used to incorporate the finite lifetime or scattering effect, which effectively replaces the delta function $\delta(\omega\pm E_{n})$ of the ${\rm Im}G$ by a Lorentzian function $1/[(\omega\pm E_{n})^{2}+\Gamma^{2}]$. In this work, the broadening parameter $\Gamma$ is assumed to have an energy dependent form, $\alpha|\omega|+3.0\times 10^{-4}$, where $\alpha=0.25$ is shown to give best agreement with experiments Choubey et al. (2017a); Wang et al. (2021); Alldredge et al. (2008). In this work, the parameters for the $t$-$t^{\prime}$-$J$ model are set to be $t^{\prime}/t=-0.3$, $J/t=0.3$. For comparison with experiments, we use $t=400$ meV, which serves as the reference for energies. The lattice sizes considered are $N_{x}\times N_{y}=48\times 48$, $48\times 96$, and $96\times 96$, which correspond approximately to magnetic fields $B\sim 11.92$, $5.96$ and $2.98$ Tesla, respectively, assuming a lattice constant $a=3.88\text{\AA}$. The number of magnetic supercells $M_{x}\times M_{y}$ are selected to ensure that $M_{x}\times N_{x}=M_{y}\times N_{y}=384$.

To obtain vortex solutions in the RMFT calculation, the various local order parameters need to be properly initialized. For example, a point singularity or $2\pi$ vorticity has to be introduced to the phase of the pairing parameters $\Delta_{ij}$. Similarly, proper modulations of the pairing amplitudes, bond variables, and hole densities are required for the initial state of the RMFT iterations. Depending on the initial conditions, our calculations uncover several vortex solutions with nearly degenerate energies, including, among others, plain $d$-wave vortex and one with intertwined PDW/CDW orders. The many vortex solutions with comparable energies obtained in our studies are reminiscent of the multitude of bulk states with intertwined orders that compete with the uniform $d$-wave SC state Tu and Lee (2016). These intriguing results from the generalized $t$-$J$ models, including both bulk and vortex solutions, highlight the complex nature of cuprate superconductors and the importance of extrinsic effects such as disorder and impurities.

In this work we consider two particular vortex solutions, which are most relevant to experiments. The first one, to be called plain vortex in the following, can be viewed as the natural topological defects of the uniform $d$-wave SC state. The corresponding pairing field is not accompanied by either PDW or CDW orders. The second solution describes a vortex structure with a checkerboard modulations in both pairing, charge-density, and magnetization. The structures of these two vortex solutions on a $N_{x}\times N_{y}=48\times 96$ square lattice are shown in Fig. 1(a) and (b). An average hole density $\delta=0.16$ and a magnetic field $B=5.96$ T is used in the KPM-RMFT calculation. Because of the periodic boundary conditions, a configuration with two vortices separated by half the linear size in the $y$ direction was used as the initial state. Large-scale calculations are necessary here to ensure that $N_{y},N_{y}\gg\xi$, where $\xi$ is the characteristic size of vortex, hence minimizes the finite size effect.

At the top of Fig. 1(a) and (b) are the density plots of the on-site hole density $\delta_{i}$ around the vortex for plain-vortex and checkerboard-halo solution, respectively. For both types of vortices, the suppression of the pairing field at the vortex core leads to an excess hole density $\delta_{\rm core}$. Importantly, the hole density $\delta_{i}$ of the second vortex solution exhibits a checkerboard halo extending over tens of lattice constants. Next we examine the configuration of the pairing field. To that end, we first define a physical pairing order parameter, which takes into account the local $d_{x^{2}-y^{2}}$ structure and the renormalization effect:

where $\mathbf{\tau}$ represents the unit vectors $\pm\hat{\mathbf{x}}$ or $\pm\hat{\mathbf{y}}$. The corresponding 3D map of the pairing amplitude $\Delta_{i}$ is shown in the middle of Fig. 1(a) and (b) for the two different vortex solutions. The point singularity of the pairing field is demonstrated in the bottom panel of Fig. 1(a) for the case of plain vortex, which shows a clear branch cut with a $2\pi$ phase jump running along the negative $x$ direction. The amplitude of the pairing field can be well approximated as $|\Delta(\mathbf{r}_{i})|\sim\Delta_{0}\,\tanh(|\mathbf{r}_{i}|/\xi_{c})$, where the vortex size is estimated to be $\xi_{c}\sim 3a$ from the fitting. Using a lattice constant $a=3.88$Å for the Cu-O plane, this corresponds to a vortex size of $2\sim 3$ nm, which agrees well with the experimental result Edkins et al. (2019). Finally, the bottom panel of Fig. 1(b) illustrates the staggered magnetization $M_{i}$ of the checkerboard-halo solution. Taking into account the renormalization effect, the staggered magnetization is given by

where $\mathbf{\tau}$ represents the unit vectors $\pm\hat{\mathbf{x}}$ and $\pm\hat{\mathbf{y}}$, and $g^{s,xy}$, $g^{s,z}$ are the Gutzwiller factors for the exchange interactions Mi . Below we present details of these vortex structures and their experimental manifestations.

By taking into account the renormalization effect due to the strong correlation in the $t$-$t’$-$J$ model, the vortex structure in a d-wave SC state is calculated for several hole densities with magnetic field in the range of 3T to 12T. Two states with almost the same energy are used to understand puzzles found by experiments. Tunneling conductance is analyzed for the plain-vortex and CB-h states to compare with the STM experiments Machida et al. (2016); Edkins et al. (2019); Gazdić et al. (2021). The absence of ZBCP in vortex core for UD and OP cuprates Machida et al. (2016); Edkins et al. (2019) but not for heavily OD samples is easily explained. The presence of PDW states in CB-h solution will suppress the conductance peak with a subgap structure. However, even if there were no PDW presence, the ZBCP is greatly reduced in conductance spectra due to the influence of $d_{x^{2}-y^{2}}$ orbital symmetry of the copper on the STM tip. Only for OD samples, the greatly reduced ZBCP is still visible as shown in the experiment of Ref. Gazdić et al. (2021). The CB-h state in the vortex halo has a bidirectional PDW state with modulation period of 8a and it shows clear $s$+$s^{\prime}$ form factors at wave vectors $Q_{x}=(1/8,0)2\pi/a$ and $2Q_{x}$ as well as peaks in y direction as reported by Ref. Edkins et al. (2019).

The presence of both QPI at low energy and CDW at higher energy of the bipartite structure reported first by Ref Machida et al. (2016) is also a distinct property of the CB-h state. In these states the interplay between QPI and CDW depend on the hole density. For UD and OP cases, the non-dispersive CDW shows strong presence in the conductance map, $g(\mathbf{k},\omega)$, for energy around or greater than half of the SC gap but it has a QPl-like energy dispersion at lower energy. The CDW signal becomes weaker as hole density increases toward the OD regime. The QPI effect is mostly contributed by the Bogoliubov quasiparticles near the nodes of SC gap; hence they are at lower energy. By contrast, the CDW is related to the PDW state at antinodes with higher energy close to pseudogap. Such a bipartite structure in the solutions without magnetic field Tu and Lee (2019) is nicely reflected in the vortex structure.

Just as reported in Ref. Edkins et al. (2019) our CB-h state has stronger $s$+$s^{\prime}$ form factors than $d$ form factors. However our state has the presence of SDW order intertwined with PDW and CDW orders while the STM experiment did not provide the magnetic information. If there were no SDW order, the $d$ form factors will dominate (This part is not included in this paper). Although there are several experiments indicating the magnetic signal inside the vortex Ref. Lake et al. (2001, 2002), more direct evidence for the presence of SDW is welcome.

Our result provided a new insight about the coupling between the bidirectional PDW states or CB with a vortex. The pairing order parameter inside a vortex is a function of the position, its DFT are FFLO orders $\langle c_{\mathbf{k}+\mathbf{Q}/2,\uparrow}c_{-\mathbf{k}+\mathbf{Q}/2,\downarrow}\rangle$. But for usual superconductors, the long coherence length will have a vortex with large radius, thus only small-$\mathbf{Q}$ FFLO orders are important. However, for cuprates, the strong correlation makes the vortex very small with a radius of three unit cells only. Thus large-$\mathbf{Q}$ FFLO orders are also important. Thus vortex could then easily induce the specific-$\mathbf{Q}$ PDW states already preferred in the cuprates without magnetic field.

In our previous works Tu and Lee (2016), we have shown that unidirectional PDW states with different periods all have similar energies in this approach using RMFT for the $t$-$J$ model. This is also found by more accurate numerical methods Zheng et al. (2017) for the Hubbard model. Thus it is difficult to identify what value of the period is more preferred for the PDW ground state from numerical calculations. Or the preference of particular period like $8a$ is just not included in the model. A recent experiment has indicated that the electron-phonon coupling may favor CDW to have a period $4a$ than $6a$ Huang et al. (2021). A consequence of the presence of FFLO orders in the vortex is that the heavily doped samples, which have no PDW states reported so far and also not favored in the theoretical calculationsTu and Lee (2016); Choubey et al. (2017a), will be induced by the magnetic field.