Superconducting Orbital Magnetoelectric Effect and its Evolution across the Superconductivity Normal Metal Phase Transition

Superconducting magnetoelectric effect, which is the current-induced magnetization in a superconductor, mainly focused on the spin magnetization in previous studies, but ignore the effect of the orbital magnetic moments carried by the paired Bloch electrons. In this work, we show that orbital magnetic moments in superconductors can induce large orbital magnetization in the presence of a current. We constructed a unified description for the current-induced spin and orbital magnetization across the superconductivity normal metal phase transition. We find that in a superconductor with uniform pairing, the current-induced magnetization at a given current density is the same as that in its normal metal state, while with the nonuniform superconducting pairing, the current-induced magnetization exhibits an abrupt change in magnitude near the superconductivity normal metal phase transition. Importantly, our theory predicts the orbital magnetoelectric effect in superconducting twisted bilayer graphene which has paired Bloch electrons with large orbital magnetic moments and negligible spin-orbit coupling. We propose that the measurement of the current-induced orbital magnetoelectric effect can be used to detect the possible nonuniform pairings in twisted bilayer graphene and other newly discovered superconductors with non-trivial Berry curvatures.

Introduction. Superconducting magnetoelectric effect is the current-induced magnetization in the superconducting state of a material. In previous studies, it mainly focused on the current-induced spin magnetization which arises from the spin-orbit coupling (SOC) in noncentrosymmetric superconductors [1, 1-6, 8, 9]. Besides the SOC, nonzero Berry curvature can also arise due to the inversion symmetry breaking [10] and has effect on the superconducting state [11]. It is known that Berry curvature introduces orbital magnetic moments to the Bloch electrons [10,[12][13][14][15][16], and the orbital magnetic moments are the source of the current-induced orbital magnetization in a normal metal [2,18]. In the superconducting state, the orbital magnetic moments from the paired Bloch electrons are involved in Cooper pairs, but the effect of orbital magnetic moments has never been studied in the superconducting magnetoelectric effect. It raises the problem of how the current-induced orbital magnetization would evolve in the phase transition from the normal metal state to the superconducting state.
In this work, we show that in superconductors where the paired Bloch electrons carry orbital magnetic moments, applying current can generate orbital magnetization. Importantly, we constructed a unified description for the current-induced spin and orbital magnetization that is applicable in both regions across the superconductor-normal metal phase transition. At the phase transition, the current-induced magnetization in the normal metal state is found to have smooth connection to that in the superconducting state when the pairing is uniform on the Fermi surfaces, while the magnetization exhibits an abrupt change in the nonuniform pairing case. The abrupt change is further ascribed to the nonuniform excitation of quasiparticles controlled by the nonuniform pairing on the Fermi surfaces.
To demonstrate the utility of our theory, we study the case of twisted bilayer graphene (TBG) and predict that the TBG exhibits current-induced orbital magnetization in the superconductivity region. The TBG has recently been observed to exhibit both superconductivity [19][20][21][22][23][24] and orbital magnetism [24][25][26][27][28]. Importantly, the superconductivity and orbital magnetism have been observed in the same TBG sample (at different filling factors) [24], indicating that the Cooper pairs are formed by electrons carrying finite orbital magnetic moments. Here we find that the orbital magnetic moments from the paired Bloch electrons can give rise to current-induced orbital magnetization in superconducting TBG. More importantly, we point out that the measurement of current-induced orbital magnetization across the superconductor-normal metal phase transition in TBG can test whether superconducting TBG has uniform or nonuniform pairing order parameters. Besides being applicable to TBG, our theory is generally applicable to a large number of noncentrosymmetric superconductors with finite Berry curvatures.
A unified description for the magnetoelectric effect. Both the normal and superconducting states of a material can be described by the Bogliubov-de Gennes Hamiltonian [29] Here H 0 (k) is the normal state Hamiltonian matrix, ∆ (k) is the pairing matrix and c † k (c k ) is the creation (annihilation) operator that includes multiple components for all the orbital and spin degrees of the system. In the normal state, a Bloch electronic state |φ ν,k with energy J (a) (b) J quasiparticle Cooper pair FIG. 1: Applying a current to a normal metal (a) and a superconductor (b). The blue and red circle correspond to the Fermi surface with and without current respectively. In a normal metal, applying a current makes more electronic states with momentum parallel to the applied current density J occupied at the Fermi energy. The occupation is schematically illustrated by the blue dot line in (a). As Bloch electrons carry magnetic moments, the redistribution at the Fermi energy gives rise to net magnetization. In a superconductor, the Bloch electrons of net momentum get paired to form Cooper pairs of net momentum. The magnetic moments carried by the paired Bloch electrons are involved in the pairing condensation and can generate a net magnetization. At finite temperature, there are quasiparticle excitations in the superconducting state, so the excited quasiparticles also contribute to the current and current-induced magnetization at 0 < T < Tc. The black arrow denotes the total magnetic moments carried by the electrons moving forward or backward, and the amplitude of the magnetic moment is schematically represented by the size of the arrow. ξ ν,k = φ ν,k | H 0 (k) |φ ν,k carries the spin magnetic moment S ν,k = φ ν,k | 1 2 µ B gσ |φ ν,k , with µ B being the Bohr magneton and g the Lande g factor. The SOC in the material is manifested by the locking of the spin magnetic moment direction with the group velocity v ν,k = ∂ k ξ ν,k . In the absence of inversion symmetry, nonzero Berry curvature Ω ν,k = i ∂ k φ ν,k | × |∂ k φ ν,k can arise and endow the Bloch electronic state with finite orbital magnetic moment m ν,k = ie [10,12,13]. Therefore, the total magnetic moment of a Bloch electron is M ν,k = S ν,k + m ν,k . In the superconducting state with time reversal symmetry, as a Bloch electronic state |φ ν,k always has a partner |φ ν,−k with the same energy ξ ν,k = ξ ν,−k , the two Bloch electronic states can pair and give rise to the pairing order parameter ∆ ν,k ∼ φ ν,−k φ ν,k [29,30], yielding the Bogliubov quasiparticle spectrum ν,k = ξ 2 ν,k + |∆ ν,k | 2 . Therefore, in a noncentrosymmetric superconductor with finite Berry curvature, the paired electrons generally carry both spin and orbital magnetic moments even though the net magnetization is zero in the absence of a current. However, as we show below, applying a current can induce net spin and orbital magnetization.
Applying a current to a metal or superconductor can be described by introducing a U (1) gauge field to the original Hamiltonian H → H + δH where with the spatial components denoted by i, j = x, y, z, and A (r, t) = q A (q, t) e iq·r being the vector gauge potential. The perturbation δH changes the distribution of Bloch electrons at the Fermi energy, and also changes the original Fermi surfaces, as is shown in Fig.  1 (a). In the normal metallic state, the redistribution of Bloch electrons results in an imbalance of the total magnetic moments of the occupied states, so a net magnetization arises by applying a current. In the superconducting state, the Bloch electrons on the Fermi surfaces with net momentum q are paired into Cooper pairs, which have the form ∆ ν,k,q ∼ φ ν,−k+ 1 2 q φ ν,k+ 1 2 q . The Cooper pairs with net momentum generate supercurrent, and also induce net magnetization that comes from the spin and orbital magnetic moments in the paired Bloch electrons. Besides the Cooper pairs, quasiparticles with net momentum are excited at finite temperature and also contribute to the net magnetization. Therefore the net magnetization in a superconductor comes from two aspects: one is the supercurrent and the other is the quasiparticle current, as is schematically illustrated in Fig. 1 (b). Importantly, the way for a supercurrent to induce net magnetization is in analogy to that in the normal metallic state: in the superconductivity the supercurrent accumulates net spin and orbital magnetic moments in the pairing condensation, while in the normal metal the current populates net magnetic moments at the Fermi energy.
To calculate the current-induced magnetization in both the normal and the superconducting states, we apply the linear response theory and obtain the current-induced bulk magnetization as [29] This magnetization is associated with the current density Here, k ≡ BZ dk/ (2π) d with d being the dimension, τ is the effective scattering time, and f ( ) is the Fermi Dirac distribution function. In the derivation of Eq. 3 and 4 we applied the Coulomb gauge so the electric field takes the form E = iωA [31]. The bulk magnetization derived in Eq. 3 involves both the spin and orbital magnetization, and it is applicable in both the normal metal and the superconductivity region. In the normal metal region with zero pairing gap ∆ ν,k = 0, Eq. 4 is the standard Drude formula to describe the current density under an applied electric field E, and the current-induced magnetization in Eq. 3 recovers the magnetoelectric susceptibility of a metal [2,18]. In the superconductivity region with finite pairing gap |∆ ν,k |, the first terms in Eq. 3 and 4 describe the net magnetization and the associated current density from the excited quasiparticles. As the quasiparticle excitations are suppressed by the pairing gap, the current and magnetization that come from the quasiparticles approach to zero at low temperature. The second terms in Eq. 3 and 4 describe the magnetization and supercurrent density from Cooper pairs. The supercurrent in Eq. 4 arises from the vector gauge field A, which twists the phase of pairing condensation. The second term in Eq. 3 is the net magnetization of pairing condensation induced by the supercurrent. Importantly, as the temperature decreases, the supercurrent and the supercurrent induced magnetization become dominant.
For the current-induced magnetization described in Eq. 3 and 4, the vector fields E and A can further be replaced by the current density J , which results in the following expression: Here the susceptibility tensor α ij describes the magnetization induced by a unit current density and can be nonzero only for crystals with point group symmetry belonging to one of the 18 gyrotropic point groups [9,32]. The tensors γ ij ,ṽ jk that are used to calculate α ij are Importantly, at zero temperature , with k F being the wave vector on the Fermi surfaces, so the induced magnetization M at a given current density J at T = 0 K is the same regardless of whether the system is superconducting or not. The susceptibility α ij at T = 0 K is only determined by the group velocity and total magnetic moments of Bloch electrons on the Fermi surfaces. It indicates that for a normal metal that has current-induced orbital magnetization, the orbital magnetization can also arise from applying supercurrent in its superconductivity region. In the temperature range 0 < T < T c (T c is the critical temperature for a superconductor), as the excited quasiparticles come into play a role in the currentinduced magnetization, the pairing gap that controls the quasiparticle excitations becomes another factor that can affect the magnetization at given current density.
Nonuniform pairing induced abrupt change in the current-induced magnetization. The effect of quasiparticle excitations on the current-induced magnetization is the most prominent at the temperature near the superconductor-normal metal phase transition, where the small pairing gap near T c allows a large number of quasiparticle excitations to coexist with Cooper pairs. From the first term in Eq. 3 and 4 it is known that the current from the excited quasiparticles and the affiliated magnetization inherit from those in the normal metal state. In the case of uniform pairing gap |∆ ν,kF | = ∆ 0 , since the suppression of quasiparticle excitations is uniform on the Fermi surfaces, the quasiparticle current and the affiliated magnetization are produced in the same ratio as the current and the current-induced magnetization in the normal metal state. In contrast, when the pairing gap ∆ ν,kF is nonuniform on Fermi surfaces, the ratio of quasiparticle current and the affiliated magnetization differs from that in the normal metal. We know from the T = 0 K result that the ratio of the supercurrent and the supercurrent-induced magnetization is the same as the ratio in the normal metal state. Therefore, across the superconductor-normal metal phase transition, the current-induced magnetization is exactly the same in the uniform pairing case, while an abrupt change in magnitude can happen in the nonuniform pairing case.
The abrupt change of current-induced magnetization at T c in the nonuniform pairing case can be deduced from the Taylor expansion of γ ij ,ṽ ij in terms of ∆ ν,k , which yields [29], so the dependence on the paring ∆ ν,kF arises abruptly in the susceptibility α ij once T < T c . When the pairing is nonuniform on the Fermi surfaces, the susceptibility α ij thus differs at the two sides of T c . Such abrupt change of current-induced magnetization across T c is an indicator of the nonuniform pairing in the superconducting state.
The orbital magnetoelectric effect in superconducting twisted bilayer graphene. In previous studies, currentinduced magnetization in the superconducting state only involves the spin magnetization caused by SOC [1, 1-6, 8, 9]. Recently, superconductivity was observed in TBG which has negligibly small SOC. Here, we predict that a large orbital magnetization can be induced by a current in superconducting TBG. More importantly, the behavior of the current-induced orbital magnetization in TBG near T c provides further information about the pairing gap in TBG.
The TBG that shows superconductivity in experiment lies on the hBN substrate, so the coupling of the bottom graphene layer with the hBN can gap out the Dirac points and generate finite Berry curvature [6,[34][35][36]. Due to the mismatch between the hBN and the bottom graphene layer, the equilateral triangle Moiré superlattice has been observed to get deformed [37][38][39], indicating that the C 3 symmetry of the TBG is generally broken by the strain from the hBN substrate. The resulting TBG system has the lowest C 1 symmetry with finite Berry curvature, so an in-plane current can give rise to an outof-plane orbital magnetization M z = α zx J x + α zy J y in the normal state [40]. Since at the T = 0 K the currentinduced magnetization in the normal state is the same as the supercurrent-induced magnetization, the currentinduced orbital magnetization in TBG is always maintained in its superconducting state. To demonstrate the in-plane current-induced out-of-plane orbital magnetization in a superconducting TBG on hBN, we construct the continuum model for a TBG with a twist angle of 1.6 • [4,5,43]. In the TBG, there is a uniaxial strain applied along the zig-zag direction in the bottom graphene layer to have it stretched by 0.1% (the detailed description for the strained TBG can be found in the Supplemental Material [29]). With the massive Dirac gap of 34meV [44,45], the Bloch electrons in the valence Moiré flat band carry orbital magnetic moments up to 80 µ B as shown in Fig. 2 (a). Around 1/2 filling in the valence Moiré flat band, the in-plane supercurrent induced The current-induced out-of-plane orbital magnetization across the superconductor-normal metal phase transition. The applied in-plane current density is 1nA/nm and the direction is along the optimal direction labelled in Fig. 2 (b). The out-of-plane orbital magnetization exhibits an abrupt jump in nonuniform pairing case while it shows smooth transition across Tc in uniform paring case.
out-of-plane orbital magnetization at T = 0K is directly computed through Eq. 6 and 7. Assuming the in-plane supercurrent density to be 1nA/nm, the out-of-plane orbital magenetization as a function of the in-plane supercurrent direction is shown in Fig. 2 (b). The largest orbital magnetization induced by the in-plane current along the optimal direction reaches the order of 10 −4 µ B /nm 2 . This magnitude of orbital magnetization is comparable to the current induced spin polarization in large Rashba SOC materials under E = 10 3 ∼ 10 4 V/m, such as Au (111) surfaces and Bi/Ag bilayers [46,47]. Near the superconductor-normal metal phase transition, the proportion of orbital magnetization that comes from the quasiparticle current becomes considerable, so the pairing gap that controls the quasiparticle excitations will affect the current-induced orbital magnetization near T c .
For the superconducting TBG on hBN, a possible nonuniform singlet pairing can be approximated as ∆ ν,k /∆ 0 = 1+λ {cos (k ·ã 1 ) + cos (k ·ã 2 ) + cos [k · (ã 1 −ã 2 )]} [48][49][50][51][52], with ∆ 0 = 1.76k B T c tanh 1.78 T c /T − 1 andã 1 , a 2 being the primitive lattice vectors of the Moiré superlattice under strain [29]. The nonuniform pairing gap at λ = 0.5 is plotted in the mini-Brillouin zone shown in Fig. 3 (a), and the orbital magnetization induced by the current along the optimal direction clearly exhibits an abrupt jump at T c as seen in Fig. 3 (b). In contrast, in the uniform pairing case which has λ = 0, the current-induced orbital magnetization has smooth connection across T c . This is consistent with our analysis that an abrupt change of the current-induced magnetization at T c indicates the pairing is nonuniform on the Fermi surfaces.
Discussion. Experimentally, the magnetization that involves orbital magnetic moments polarization in the normal state of materials has been observed through the optical Kerr effect [53], nuclear magnetic resonance measurements [54], and the superconducting quantum interference device (SQUID) [28], and these techniques can all be used to detect the magnetization in the superconducting state [55][56][57]. In the recent direct image of orbital magnetism in TBG [28], it can be deduced from the measurement that the surface magnetism at the order of 10 −4 µ B /nm 2 corresponds to the generated static magnetic field around 1nT, so a resolution of 0.1nT SQUID device will be able to detect the orbital magnetization induced by the current density at 1nA/nm in superconducting TBG. The current-induced orbital magnetization and the possible abrupt change at T c is proportional to the applied current density in the material, so the critical current density in the superconducting state sets the upper limit of the current-induced orbital magnetization below T c . As a result, superconductors with larger pairing gap that have larger critical current density would be preferred to make the effect more measurable.
It is important to note that our theory applies to a large number of recently discovered superconductors with Berry curvatures and negligible SOC, such as a trilayer graphene Moiré superlattice on hBN [58][59][60][61], bilayer graphene/hBN superlattices [62], and twisted double bilayer graphene [63,64]. Moreover, our theory also applies to non-centrosymmetric superconductors with strong SOC and finite Berry curvatures such as chiral crystals,

BOGLIUBOV-DE GENNES HAMILTONIAN AND THE GREEN'S FUNCTION
The generic Hamiltonian that can describe the normal state of a system takes the form with c † ν,k (c ν,k ) being the creation (or annihilation) operator, H 0,νν (k) being the element of the Hamiltonian matrix H 0 (k), and ν = 1, 2, 3, . . . denoting the spin and orbital index of the system. The Bogliubov-de Gennes Hamiltonian that can describe the superconductivity state of a system is where . . , and∆ (k) is the pairing matrix. We know that the normal state Hamiltonian H 0 (k) can be diagonalized by an unitary transformation U (k) as with U † (k) U (k) = 1 and the corresponding eigen states φ ν,k = U νν (k) c ν ,k , φ † ν,k = c † ν ,k U * ν ν (k). Applying the unitary transformation to the Bogliubov-de Gennes Hamiltonian then yields Importantly, after the unitary transformation, we obtain In such eigen-band basis φ † 1,k , φ † 2,k , ..., φ † ν,k , φ 1,−k , φ 2,−k , ..., φ ν,−k , the intraband pairing φ ν,−k φ ν,k gives an energetically favourable pairing phase, while other interband pairing will require extra attractive interaction to overcome the energy (momentum) difference [S1]. In the energetically favourable pairing phase, the pairing matrix respects H 0 (k)∆ k −∆ k H * 0 (−k) = 0 [S1], and the pairing matrix is therefore diagonalized simultaneously by the unitary transformation as As a result, the Bogliubov-de Gennes Hamiltonian in the eigen-band basis takes the form We define the Green's function for the Bogliubov-de Gennes Hamiltonian as , and G 0 (k, iω n ) = G e (k, iω n ) F (k, iω n ) and then we apply the unitary transformation to the Green's function: which gives (S13)

LINEAR RESPONSE THEORY FOR THE CURRENT-INDUCED MAGNETIZATION
Applying current introduces a perturbation to the Bogliubov-de Gennes Hamiltonian, so the free-energy density of the system gets changed. The change of free-energy density can be obtained through expanding the electromagnetic gauge fields, and it can be expressed as Here, the terms T s give the spin magnetoelectric susceptibility in both the superconductivity and normal metal regions, while the polarization tensor Π ij (q, iω m ), where is responsible for the current and current-induced orbital magnetization in both the superconductivity and normal metal regions.

Polarization Tensor, Current, and the Current-induced Orbital Magnetization
Substituting Eq. S10, S11, S12, S13 into Eq. S18, we can obtain the polarization tensor Π ij (q, iω m ), where with the static part Π sta ij (q) and the dynamic part Π dyn ij (q, iω m ), where The static and dynamic part of the polarization tensor can then be expanded in terms of q to the linear order as with The intra-band contribution to Π sta,(0) ij is Π sta,(0),intra ij For the inter-band contribution, we assume that the pairing gap |∆ ν,k | is much smaller than the inversion symmetry breaking-induced band splitting, so the inter-band term can be approximated by that used in the normal metal phase with ∆ ν,k = 0: We also know from the gauge invariance that with the group velocity v ν,k = ∇ k ξ ν,k = φ ν,k | ∇ k H 0 (k) |φ ν,k . For the term Π dyn,(0) ij that is iω m dependent, we consider the dominant intra-band contribution and take the analytic continuation iω m → ω + i τ −1 to get The intra-band contribution to Q sta ijl [S2, S3] is with m ν,k = ie 2 ∂ k φ ν,k |×[H 0 (k) − ξ ν,k ] |∂ k φ ν,k being the orbital magnetic moment. For the inter-band contribution, given that the pairing gap |∆ ν,k | is much smaller than the band splitting, it can also be approximated by that used in the normal metal phase in a way that is similar to Eq. S23. We know that the summation of the intra-band and inter-band static terms in the normal metal phase gives Q sta ijl = e 2 2 BZ dk (2π) d ν 2f ( ν,k ) ∇ k · ( ν,k Ω ν,k ), which is confirmed to be zero [S2, S3]. As a result, the inter-band contribution to Q sta ijl is For the dynamic part Q dyn ijl , we consider the dominant intra-band term and use the analytical continuation iω m → ω + i τ −1 as well, and then we can obtain Eventually, the polarization tensor Π ij (q, ω) can be approximated: From Eq. S15, we can then obtain the current density: Interestingly, we notice that the free energy density related to Q sta ijl and Q dyn ijl can be further simplified as Then, we derive the current-induced orbital magnetization: Here the electric field has been expressed in terms of the vector gauge potential as E = iωA.

Spin Magnetoelectric Susceptibility Calculation
The spin magnetoelectric susceptibility T s ij is calculated as: We further sum over the Matsubara frequency and decompose the spin magnetoelectric susceptibility into the static (iω m = 0) and dynamic (iω m = 0) parts as with T s,sta and The intra-band contribution to T s,sta ij is T s,sta,intra For the inter-band contribution, similarly, we assume that the pairing gap |∆ ν,k | is much smaller than the band splitting so that it can be approximated by that in the normal metal phase. We also know that the summation of the intra-band and inter-band terms in the static limit will vanish, so the inter-band contribution can be written as We know that for a Bloch state, the spin magnetic moment is S ν,k = φ ν,k | 1 2 µ B gσ |φ ν,k , so the static part of the spin magnetoelectric susceptibility can be calculated as can be obtained similarly following the above procedure: For the dynamic part that is iω m dependent, we consider the dominant intra-band contribution and obtain T s,dyn ij (iω m ) can be obtained similarly following the above procedure as well: Eventually, we get the spin magnetoelectric susceptibility as We consider 42 sites in the hexagonal reciprocal lattice so that h η (k) is an 84×84 matrix in the calculation for the energy dispersion ξ ν,s,η,k . For a specific band with index ν, the Bogliubov de-Gennes Hamiltonian in the eigen-band basis φ † ν,↑,η,k , φ † ν,↓,η,k , φ ν,↓,−η,−k , φ ν,↑,−η,−k can be written as with the singlet pairing ∆ k = ∆ 0 + λ∆ 0 {cos (k ·ã 1 ) + cos (k ·ã 2 ) + cos [k · (ã 1 −ã 2 )]}.