text:
# Preserving the gauge invariance of meson production currents in the presence of explicit final-state interactions

## I Introduction

Gauge invariance is one of the central issues when attempting to describe how photons interact with hadronic systems. Concentrating on the simplest case—pion photoproduction with real or virtual photons, gauge invariance can easily be shown to follow if the $`\pi N`$ and $`\gamma N`$ problems are treated completely and consistently on an equal footing . In practice, however, one often needs to revert to some approximate treatment of one or more of the contributing reaction mechanisms and this usually leads to a violation of gauge invariance. To restore it, the neglected reaction mechanisms must be approximated by auxiliary currents constructed such that the gauge-invariance-violating contributions to the four-divergence of the total production amplitude are cancelled. Such a procedure cannot be unique, of course, since one may always add arbitrary transverse currents without affecting the four-divergence.

At the tree-level, where one does not resolve the internal mechanisms entering the interaction current (cf. Fig. 1), various recipes exist to preserve gauge invariance. The simplest case concerns the choice of bare vertices with pseudovector coupling for the $`\pi NN`$ vertex, where the corresponding Kroll–Ruderman contact current  follows from the minimal substitution procedure. The case of extended nucleons, whose internal structure is described in terms of (phenomenological) form factors, is treated in Refs. .

In the present work, we want to go beyond the tree level and investigate how one can preserve gauge invariance if the internal structure of the interaction current is taken into account explicitly. The reaction mechanisms that enter the interaction current are summarized within the dashed box in Fig. 1. Specifically, we are interested in preserving gauge invariance in the explicit presence of hadronic final-state interactions. This is achieved by introducing auxiliary currents which cancel the gauge-invariance-violating contributions. In particular, we show how one may exploit the constraints following from the generalized Ward–Takahashi identities to construct these currents.

The present discussion is restricted to nucleons and pions only to facilitate the presentation. We do not include here possible resonances or other transition mechanisms since their couplings to the electromagnetic field are transverse and have no bearing on the question of gauge invariance. None of these restrictions are essential, however, and one may easily adapt the present formalism to accommodate more complex situations.

## II Gauge Invariance

The pion photoproduction current of the nucleon is shown in Fig. 1 . According to the diagrams in the first line of this figure, the total current $`M^\mu `$ may be broken up into four main contributions: The three Born terms due to the $`s`$-, $`u`$-, and $`t`$-channel currents stemming from the photon coupling to the three external legs of the $`\pi NN`$ vertex, and the interaction current $`M_{\mathrm{int}}^\mu `$ where the photon attaches itself to an internal leg of the $`\pi NN`$ vertex, i.e.,

$$M^\mu =M_s^\mu +M_u^\mu +M_t^\mu +M_{\mathrm{int}}^\mu .$$
(1)
While the first three contributions are relatively straightforward, the last one—as it is shown in the last two lines of Fig. 1—explicitly involves the full complexity of the internal reaction dynamics of the underlying $`\pi N`$ scattering problem summarized in Fig. 2.

For gauge invariance of the total current $`M^\mu `$ to hold true, its four-divergence must satisfy a generalized Ward–Takahashi identity

$`k_\mu M^\mu `$ $`=`$ $`[F_s\tau ]S_{p+k}Q_iS_p^1+S_p^{}^1Q_fS_{p^{}k}[F_u\tau ]`$ (3)
$`+\mathrm{\Delta }_q^1Q_\pi \mathrm{\Delta }_{qk}[F_t\tau ],`$
where $`p`$ and $`k`$ are the four-momenta of the incoming nucleon and photon, respectively, and $`p^{}`$ and $`q`$ are the four-momenta of the outgoing nucleon and pion, respectively, related by momentum conservation $`p^{}+q=p+k`$. $`S`$ and $`\mathrm{\Delta }`$ are the propagators of the nucleons and pions, respectively, with their subscripts denoting the available four-momentum for the corresponding hadron; $`Q_i`$, $`Q_f`$, and $`Q_\pi `$ are the initial and final nucleon and the pion charge operators, respectively. $`[F_x\tau ]`$ denotes the $`\pi NN`$ vertex (including coupling and isospin operators), with the subscript $`x`$ labeling the kinematic situation appropriate for the $`s`$-, $`u`$-, or $`t`$-channel diagrams appearing in Fig. 1. The vertex isospin $`\tau `$ does not commute with the charge operators and $`\tau Q_iQ_f\tau Q_\pi \tau =0`$ describes charge conservation at the vertex in a symbolic manner.

Equation (3) is easily obtained from Eq. (1) upon using the Ward–Takahashi identities  for the nucleon and pion currents,

$`k_\mu J_\mathrm{n}^\mu `$ $`=`$ $`S_{p_\mathrm{n}+k}^1Q_\mathrm{n}Q_\mathrm{n}S_{p_\mathrm{n}}^1,`$ (5)
$`k_\mu J_\pi ^\mu `$ $`=`$ $`\mathrm{\Delta }_{q_\pi +k}^1Q_\pi Q_\pi \mathrm{\Delta }_{q_\pi }^1,`$ (6)
where $`p_\mathrm{n}`$ and $`q_\pi `$ are the respective initial hadron momenta of the electromagnetic vertices, and using the fact that for gauge invariance to be true the interaction current must obey

$$k_\mu M_{\mathrm{int}}^\mu =[F_s\tau ]Q_i+Q_f[F_u\tau ]+Q_\pi [F_t\tau ].$$
(7)
### A Preserving gauge invariance

As alluded to above, the preceding relations can easily be shown to be true if the $`\pi N`$ and $`\gamma N`$ problems are treated consistently on an equal footing . In practical applications, however, approximations are inevitable which usually violate the gauge invariance.

To see how one may preserve gauge invariance in such a situation, let us explicitly write the four-divergence of the interaction current using the relevant parts of Fig. 1 as guidance. One finds

$`k_\mu M_{\mathrm{int}}^\mu `$ $`=`$ $`k_\mu m_{\mathrm{bare}}^\mu +k_\mu U^\mu G_0[F\tau ]`$ (10)
$`+XG_0\{k_\mu J_\mathrm{n}^\mu S_{p_\mathrm{n}k}[F_u\tau ]+k_\mu J_\pi ^\mu \mathrm{\Delta }_{q_\pi k}[F_t\tau ]`$
$`+k_\mu m_{\mathrm{bare}}^\mu +k_\mu U^\mu G_0[F\tau ]\},`$
where $`G_0=S_{p_\mathrm{n}}\mathrm{\Delta }_{q_\pi }`$ denotes the product of the intermediate nucleon and pion propagators, with respective four-momenta $`p_\mathrm{n}`$ and $`q_\pi `$ denoting the integration variables, and $`X`$ is the nonpolar $`\pi N`$ amplitude (see Fig. 2) which mediates the hadronic final-state interaction; $`U^\mu `$ subsumes all exchange currents and $`m_{\mathrm{bare}}^\mu `$ is the bare contact current.

For gauge invariance to hold true, the bare current must satisfy the condition

$`k_\mu m_{\mathrm{bare}}^\mu =[f_s\tau ]Q_i+Q_f[f_u\tau ]+Q_\pi [f_t\tau ],`$ (11)
where $`[f_x\tau ]`$ denotes the bare $`\pi NN`$ vertex the same way the notation $`[F_x\tau ]`$ was used above for the dressed vertex. This is the analog of Eq. (7) for the bare current; it is usually satisfied as a matter of course. One of the simplest nontrivial examples is the case of pure pseudovector coupling without form factors, where $`m_{\mathrm{bare}}^\mu `$ is the Kroll–Ruderman contact current ; see Eq. (27).

Combining now Eq. (10) with the necessary condition (7), and making use of Eq. (11) in the Born terms, produces

$`0`$ $`=`$ $`([F_s\tau ][f_s\tau ])Q_iQ_f([F_u\tau ]f_u\tau ])`$ (15)
$`Q_\pi \left([F_t\tau ][f_t\tau ]\right)+k_\mu U^\mu G_0[F\tau ]`$
$`+XG_0\{k_\mu J_\mathrm{n}^\mu S_{p_\mathrm{n}k}[F_u\tau ]+k_\mu J_\pi ^\mu \mathrm{\Delta }_{q_\pi k}[F_t\tau ]`$
$`+k_\mu m_{\mathrm{bare}}^\mu +k_\mu U^\mu G_0[F\tau ]\}`$
as a necessary off-shell requirement for gauge invariance to be satisfied.

In other words, as long as the basic Ward–Takahashi identities (2) for the hadron currents are true, any approximation of the full reaction mechanisms constructed in such a manner that the condition (15) is satisfied will also preserve gauge invariance as a matter of course.

In view of the arbitrariness of transverse contributions, there are of course infinitely many ways this can be achieved. The prescription we give in the following applies to the simplifying assumption that one completely omits the explicit treatment of exchange currents $`U^\mu `$. However, even if they are taken into account in some partial manner, it is a straightforward exercise to adapt the following formulation to accommodate such situations.

Omitting explicit exchange currents, we maintain gauge invariance by constructing auxiliary currents $`j_0^\mu `$ and $`j_1^\mu `$ which provide the same effect as the exchange currents $`U^\mu `$ as far as the preservation of gauge invariance is concerned.

To this end, we make the replacements

$`U^\mu G_0[F\tau ]`$ $``$ $`j_0^\mu +\mathrm{\Delta }j_0^\mu ,`$ (17)
$`XG_0U^\mu G_0[F\tau ]`$ $``$ $`j_1^\mu +\mathrm{\Delta }j_1^\mu ,`$ (18)
and demand that $`j_0^\mu `$ and $`j_1^\mu `$ satisfy

$`k_\mu j_0^\mu `$ $`=`$ $`\left([F_s\tau ][f_s\tau ]\right)Q_i+Q_f\left([F_u\tau ][f_u\tau ]\right)`$ (20)
$`+Q_\pi \left([F_t\tau ][f_t\tau ]\right)`$
and

$`k_\mu j_1^\mu `$ $`=`$ $`XG_0\{k_\mu J_\mathrm{n}^\mu S_{p_\mathrm{n}k}[F_u\tau ]+k_\mu J_\pi ^\mu \mathrm{\Delta }_{q_\pi k}[F_t\tau ]`$ (22)
$`+k_\mu m_{\mathrm{bare}}^\mu \}.`$
Clearly, if these conditions are met, then

$$k_\mu \left(\mathrm{\Delta }j_0^\mu +\mathrm{\Delta }j_1^\mu \right)=0,$$
(23)
i.e., $`\mathrm{\Delta }j_0^\mu +\mathrm{\Delta }j_1^\mu `$ is purely transverse. Without explicit treatment of exchange currents, these contributions are inaccessible and will be dropped. The resulting production current,

$`M^\mu `$ $`=`$ $`M_s^\mu +M_u^\mu +M_t^\mu +m_{\mathrm{bare}}^\mu +j_0^\mu `$ (25)
$`+XG_0\left\{M_u^\mu +M_t^\mu +m_{\mathrm{bare}}^\mu \right\}+j_1^\mu ,`$
will then satisfy the generalized Ward–Takahashi identity (3) and therefore it will be gauge invariant.

To be more specific as to how to implement the conditions (20) and (22), allowing for a mixture of pseudoscalar and pseudovector couplings, let us write the dressed $`\pi NN`$ vertex as

$$F=g_{\mathrm{ps}}\gamma _5G_{\mathrm{ps}}+g_{\mathrm{pv}}\frac{\gamma _5q/_\pi }{2m}G_{\mathrm{pv}},$$
(26)
where the indices ps and pv stand for pseudoscalar and pseudovector contributions, respectively; $`G_{\mathrm{ps}}`$ and $`G_{\mathrm{pv}}`$ denote the corresponding normalized form factors (with their strength parameters $`g_{\mathrm{ps}}`$ and $`g_{\mathrm{pv}}`$ adding up to the physical coupling constant, $`g_{\pi NN}=g_{\mathrm{ps}}+g_{\mathrm{pv}}`$), $`q_\pi `$ is the four-momentum of the pion, and $`m`$ the nucleon mass. The bare vertex $`f`$ is given by the same equation with all $`G`$’s removed, and the corresponding bare current,

$$m_{\mathrm{bare}}^\mu =g_{\mathrm{pv}}\frac{\gamma _5\gamma ^\mu }{2m}Q_\pi \tau ,$$
(27)
is just the usual Kroll–Ruderman contact term .

Equation (20) now reads explicitly

$`k_\mu j_0^\mu `$ $`=`$ $`\gamma _5\left(G_sg\right)\tau Q_i`$ (32)
$`+\gamma _5\left(G_ug\right)Q_f\tau `$
$`+\gamma _5\left(G_tg\right)Q_\pi \tau `$
$`g_{\mathrm{pv}}\gamma _5{\displaystyle \frac{k/}{2m}}\left(G_{\mathrm{pv},s}\tau Q_iG_{\mathrm{pv},u}Q_f\tau \right)`$
$`k_\mu m_{\mathrm{bare}}^\mu ,`$
where

$$G_x=g_{\mathrm{ps}}G_{\mathrm{ps},x}+g_{\mathrm{pv}}\frac{p/p/^{}}{2m}G_{\mathrm{pv},x}$$
(33)
(with $`x=s`$, $`u`$, $`t`$) denotes the kinematic situations in which the vertex functions $`G_{\mathrm{ps}}`$ and $`G_{\mathrm{pv}}`$ appear; $`g`$ is given by the same equation with all $`G`$’s removed.

Note that all the terms containing $`g`$ in Eq. (32) add up to zero. We may therefore replace $`g`$ by an arbitrary function $`\widehat{F}`$ with impunity. Furthermore, using the Mandelstam variables $`s=(p+k)^2`$, $`u=(p^{}k)^2`$, and $`t=(qk)^2`$, we may then rewrite the resulting equation as

$`k_\mu j_0^\mu `$ $`=`$ $`k_\mu \{{\displaystyle \frac{(2p+k)^\mu }{sp^2}}\gamma _5(G_s\widehat{F})\tau Q_i`$ (38)
$`{\displaystyle \frac{(2p^{}k)^\mu }{up^2}}\gamma _5\left(G_u\widehat{F}\right)Q_f\tau `$
$`{\displaystyle \frac{(2qk)^\mu }{tq^2}}\gamma _5\left(G_t\widehat{F}\right)Q_\pi \tau `$
$`g_{\mathrm{pv}}{\displaystyle \frac{\gamma _5\gamma ^\mu }{2m}}\left(G_{\mathrm{pv},s}\tau Q_iG_{\mathrm{pv},u}Q_f\tau \right)`$
$`m_{\mathrm{bare}}^\mu \},`$
which allows us to put

$`j_0^\mu `$ $`=`$ $`{\displaystyle \frac{(2p+k)^\mu }{sp^2}}\gamma _5\left(G_s\widehat{F}\right)\tau Q_i`$ (42)
$`{\displaystyle \frac{(2p^{}k)^\mu }{up^2}}\gamma _5\left(G_u\widehat{F}\right)Q_f\tau `$
$`{\displaystyle \frac{(2qk)^\mu }{tq^2}}\gamma _5\left(G_t\widehat{F}\right)Q_\pi \tau `$
$`g_{\mathrm{pv}}{\displaystyle \frac{\gamma _5\gamma ^\mu }{2m}}\left(G_{\mathrm{pv},s}\tau Q_iG_{\mathrm{pv},u}Q_f\tau \right)m_{\mathrm{bare}}^\mu .`$
Comparing with Eq. (25), note that the last two terms of this gauge-invariance-preserving current cancel the bare Kroll–Ruderman term and replace it by a dressed one, where instead of the pion charge $`Q_\pi \tau `$, there is now a dressing term $`G_{\mathrm{pv},s}\tau Q_iG_{\mathrm{pv},u}Q_f\tau `$. This dressing term expresses the pion charge in terms of the nucleon charges modified by hadronic form factors, and, in general, it will be non-zero even if the pion is uncharged.

We emphasize that the transition from Eq. (32) to (42) is not unique, of course, since one may add a divergence-free current to $`j_0^\mu `$ without changing the necessary condition (32). In fact, the replacement of $`g`$ by $`\widehat{F}`$ amounts to the addition of such a transverse current (which in turn may be understood as a phenomenological way of getting a handle on the neglected transverse current $`\mathrm{\Delta }j_0^\mu +\mathrm{\Delta }j_1^\mu `$). At this stage, $`\widehat{F}`$ is completely undetermined. Below, when considering the relationship of the present results to existing tree-level approaches, we will discuss some specific choices.

We also note that the terms appearing in Eq. (42) do not introduce any new singularities into the amplitude. This is easily illustrated for the example of the $`s`$-channel pole diagram $`M_s^\mu `$, viz.

$`M_s^\mu `$ $`=`$ $`[F_s\tau ]{\displaystyle \frac{p/+k/+m}{sm^2}}\gamma ^\mu Q_i+M_{\mathrm{t},s}^\mu `$ (43)
$`=`$ $`\gamma _5G_s{\displaystyle \frac{(2p+k)^\mu }{sm^2}}\tau Q_i+g_{\mathrm{pv}}{\displaystyle \frac{\gamma _5k/}{2m}}G_{\mathrm{pv},s}{\displaystyle \frac{(2p+k)^\mu }{sm^2}}\tau Q_i`$ (45)
$`+\gamma _5G_{s,q/}{\displaystyle \frac{k^\mu \gamma ^\mu k/}{sm^2}}\tau Q_i+M_{\mathrm{t},s}^\mu ,`$
where $`G_{s,q/}`$ is given by Eq. (33) with $`p/p/^{}`$ replaced by $`q/`$; $`M_{\mathrm{t},s}^\mu `$ splits off the transverse part of the electromagnetic nucleon current given in Eq. (50) below. The decomposition given here makes it immediately obvious that the effect of adding the $`s`$-channel term, with $`p^2=m^2`$, of Eq. (42) to this expression is to replace $`G_s`$ in the first term here by $`\widehat{F}`$. The same happens also for the $`u`$\- and $`t`$-channel diagrams. In other words, apart from providing a dressed Kroll–Ruderman term, the effect of the gauge-invariance preserving current (42) is to provide a common form factor $`\widehat{F}`$ for some (but not all) of the contributions originally containing individual form factors $`G_s`$, $`G_u`$, and $`G_t`$, without changing the original singularities of the Born diagrams.

Next, to construct the current $`j_1^\mu `$, we recall that the gauge-invariant nucleon and pion currents appearing in the $`u`$\- and $`t`$-channel terms of the final-state interaction contribution are given by

$`J_\mathrm{n}^\mu `$ $`=`$ $`\gamma ^\mu Q_\mathrm{n}+J_{\mathrm{t},\mathrm{n}}^\mu ,`$ (47)
$`J_\pi ^\mu `$ $`=`$ $`(2q_\pi k)^\mu Q_\pi +J_{\mathrm{t},\pi }^\mu ,`$ (48)
respectively, with transverse pieces which follow from demanding the validity of the Ward–Takahashi identities (2), i.e.,

$`J_{\mathrm{t},\mathrm{n}}^\mu `$ $`=`$ $`\left(\gamma ^\mu k^\mu {\displaystyle \frac{k/}{k^2}}\right)Q_\mathrm{n}(F_11)+i{\displaystyle \frac{\sigma ^{\mu \nu }k_\nu }{2m}}\kappa F_2,`$ (50)
$`J_{\mathrm{t},\pi }^\mu `$ $`=`$ $`\left[(2q_\pi k)^\mu k^\mu {\displaystyle \frac{k(2q_\pi k)}{k^2}}\right]Q_\pi (F_\pi 1),`$ (51)
where $`F_1`$ and $`F_2`$ respectively are the electromagnetic Dirac and Pauli form factors of the nucleon, with $`\kappa `$ being its anomalous magnetic moment, and $`F_\pi `$ is the electromagnetic form factor of the pion. Their four-divergence is given by

$`k_\mu J_\mathrm{n}^\mu `$ $`=`$ $`k_\mu \gamma ^\mu Q_\mathrm{n},`$ (54)
$`k_\mu J_\pi ^\mu `$ $`=`$ $`k_\mu (2q_\pi k)^\mu Q_\pi .`$ (55)
Inserting this into Eq. (22), we may then extract the gauge-invariance-preserving current as

$`j_1^\mu `$ $`=`$ $`XG_0\{\gamma ^\mu Q_\mathrm{n}S_{p_\mathrm{n}k}[F_u\tau ]`$ (57)
$`+(2q_\pi k)^\mu Q_\pi \mathrm{\Delta }_{q_\pi k}[F_t\tau ]+m_{\mathrm{bare}}^\mu \}.`$
Its effect on the final-state interaction part of the production amplitude (25) is seen to simply cancel the bare current and to reduce the $`u`$\- and $`t`$-channel contributions to their respective transverse pieces, i.e.,

$`XG_0\left(M_{\mathrm{t},u}^\mu +M_{\mathrm{t},t}^\mu \right)`$ $`=`$ $`XG_0\left\{M_u^\mu +M_t^\mu +m_{\mathrm{bare}}^\mu \right\}+j_1^\mu `$ (58)
$`=`$ $`XG_0\{J_{\mathrm{t},\mathrm{n}}^\mu S_{p_\mathrm{n}k}[F_u\tau ]`$ (60)
$`+J_{\mathrm{t},\pi }^\mu \mathrm{\Delta }_{q_\pi k}[F_t\tau ]\}.`$
The entire contribution from the final-state interaction, therefore, is purely transverse.

### B Relation to Ohta’s and Haberzettl’s tree-level prescriptions

To make the connection with tree-level approaches, we need to switch off all final-state interactions and put $`X=0`$ in Eqs. (22) and (25). The conditions to be satisfied, therefore, are Eq. (20) and Eq. (22) in the form

$$k_\mu j_1^\mu =0.$$
(61)
In other words, simply putting $`j_1^\mu =0`$ satisfies all gauge-invariance constraints at the tree level.

Both Ohta’s  and Haberzettl’s  prescriptions for preserving gauge invariance can be understood as different choices for the function $`\widehat{F}`$ in Eq. (42).

Ohta’s approach, based on a particular application of the minimal substitution procedure, finds

$$\widehat{F}=g_{\mathrm{ps}}G_{\mathrm{ps}}(q,p^{},p)+g_{\mathrm{pv}}\frac{p/p/^{}}{2m}G_{\mathrm{pv}}(q,p^{},p),$$
(62)
where the external hadron momenta of the photoproduction current—a four-point function—appear here as the momenta of the $`\pi NN`$ vertex—a three-point function. Since the momenta satisfy $`q+p^{}=p+k`$, this mismatch corresponds to an unphysical region of the $`\pi NN`$ vertex (which leads to the problems discussed in Ref. ). Only in the infrared limit of $`k=0`$, this mismatch is resolved and then this choice prevents the current (42) from being singular at $`k=0`$.

In Haberzettl’s prescription, the function $`\widehat{F}`$ is a linear combination of the three kinematical situations in which the $`\pi NN`$ vertices appear in the Born terms, i.e., \[cf. Eq. (33)\]

$$\widehat{F}=a_sG_s+a_uG_u+a_tG_t,$$
(63)
with coefficients constrained by $`a_s+a_u+a_t=1`$, which may be fixed according to prejudice or used as free fit parameters. In contrast to Ohta’s choice, this does not require any unphysical values for the $`\pi NN`$ form factors in practical applications, and it also has a well-behaved infrared limit. In direct comparisons to Ohta’s, this prescription is found to provide better agreement with the experimental data  .

## III Summary and Discussion

In summary, we have treated here the electromagnetic production current for mesons off the nucleon for both real and virtual photons. We used the constraints following from requiring the validity of the generalized Ward–Takahashi identities to construct auxiliary current pieces that ensure that gauge invariance is preserved even if—as it is invariably the case in practical applications—one does not treat the problem completely and consistently. The result for the production current $`M^\mu `$ obtained here in Eqs. (25), (42), and (57) can be summarized by

$`M^\mu `$ $`=`$ $`M_s^\mu +M_u^\mu +M_t^\mu +j_{\mathrm{gip}}^\mu `$ (65)
$`+XG_0\left(M_{\mathrm{t},u}^\mu +M_{\mathrm{t},t}^\mu \right),`$
where the gauge-invariance-preserving current,

$$j_{\mathrm{gip}}^\mu =j_0^\mu +m_{\mathrm{bare}}^\mu ,$$
(66)
is given via $`j_0^\mu `$ of Eq. (42).

As far as the choice of the function $`\widehat{F}`$ appearing in $`j_0^\mu `$ is concerned, the generally better results obtained with Haberzettl’s prescription at the tree-level  seem to favor the form given in Eq. (63). We emphasize, however, that other functions are possible and that, in general, as far as gauge invariance is concerned, any current $`j_0^\mu `$ that satisfies the necessary condition (32) is permitted here.

The part describing the hadronic final-state interaction due to $`X`$ does no longer contain the Kroll–Ruderman contact term since it is canceled by the current of Eq. (57). Moreover, what remains is seen to be entirely transverse since the corresponding $`u`$\- and $`t`$-channel contributions, $`M_{\mathrm{t},u}^\mu `$ and $`M_{\mathrm{t},t}^\mu `$, respectively contain only the transverse electromagnetic nucleon and pion operators of Eq. (II A). For real photons, in particular, both $`F_\pi 1`$ and $`F_11`$ vanish and therefore $`M_{\mathrm{t},t}^\mu =0`$, i.e., the $`t`$-channel does not contribute at all, and the $`u`$-channel current $`M_{\mathrm{t},u}^\mu `$ is reduced to the magnetic $`\sigma ^{\mu \nu }k_\nu `$ term from Eq. (50).

We emphasize that even though the gauge-invariant production current of Eq. (65) is an approximation to the full dynamics of the problem as summarized in Figs. 1 and 2, it is complete as far as the longitudinal components of the current are concerned. Any additional pieces must be transverse.

Let us repeat once again that the present work was restricted to pions and nucleons merely to simplify the presentation. The concepts developed here are quite general, however. An extension to other baryons and mesons, therefore, is straightforward and easily done along the lines given here. In particular, the fact that the final-state interaction contributions are purely transverse remains true if one takes into account additional intermediate hadrons ($`\mathrm{\Delta }`$, $`\rho `$, etc.) since, just as was demonstrated here for the nucleons and pions, only the transverse parts of their respective current operators survive in the final-state interaction terms.

Finally, let us point out that the present results remain equally valid whether the form factors $`F`$ and the final-state amplitude $`X`$ are obtained via some sophisticated Bethe–Salpeter-type formalism or are based on a simple phenomenological model ansatz. How these elements are obtained does not enter any of the present considerations and therefore has no bearing on the question of gauge invariance.

###### Acknowledgements.

The author gratefully acknowledges discussions with S. Krewald and K. Nakayama which precipitated the present work. This work was supported in part by Grant No. DE-FG02-95ER-40907 of the U.S. Department of Energy.


format: text

text:
# Vortex in a 𝑑-wave Superconductor at Low Temperatures

## I INTRODUCTION

Vortices in classic superconductors involve a winding of $`2\pi `$ of the order parameter phase $`\varphi (𝐫)`$ around the vortex line. Since the order parameter $`\mathrm{\Delta }_𝐤=|\mathrm{\Delta }_𝐤|e^{i\varphi }`$ in a $`d_{x^2y^2}`$-wave superconductor like the high-$`T_c`$ cuprates is also simply a complex scalar with single global phase $`\varphi `$, it was initially expected that the vortex state in the cuprates might be structurally quite similar to the textbook case. One remarkable difference was pointed out by Joynt: in the $`d`$-wave case, if subdominant pair potential components of different symmetry exist, corresponding order parameter components can be induced at $`T_c`$ by any probe which couples to gradients of the order parameter, even if the subdominant zero field “bare” critical temperatures are very small or zero. Thus the structure of a $`d`$-wave vortex generically involves admixtures of different symmetry order parameters. Secondly, as noted by Volovik, the traditional roles played by extended and localized quasiparticle states in classic superconductors are reversed in the $`d`$-wave case. In classic superconductors at low temperatures, Caroli-de Gennes-Matricon bound states in the vortex core  dominate the electronic density of states because extended states are fully gapped and therefore unoccupied. In the $`d`$-wave case, the existence of order parameter nodes inhibits the formation of the bound states (their very existence is questionable ), and populates the extended ones, which are found to dominate thermodynamics at low temperatures and magnetic fields.

Several authors  have attacked the $`d`$-wave mixed state structure problem in recent years, armed with these ideas. Early studies focussed on an isolated $`d`$-wave vortex, allowing for an induced $`s`$-wave order parameter component and solving the Bogoliubov-de Gennes equations on a lattice. Within the Ginzburg-Landau (GL) theory, similar results were obtained. The $`s`$ component was shown to have opposite winding number to that of the parent $`d_{x^2y^2}`$ in the core regions. Far away from the core center, it decays as $`1/r^2`$, and its winding number becomes 3,, implying that there are four extra vortices in the $`s`$-field at large distances from the main core.  These results were confirmed by numerical solutions of the Eilenberger equations by Ichioka et al. The possibility of an induced $`d_{xy}`$ component was also allowed for in Ref. , which concluded that the structure in this case was similar to that of the (induced) $`s`$-wave case, except that the induced order parameter at large distances was found to decay more rapidly, roughly as $`1/r^4`$, and have opposite winding number 5. Koyama and Tachiki pointed out, however, that if the calculation is done in a gauge invariant manner there is an additional term in the free energy not found by Ichioka et al., involving a Zeeman coupling of the field to an intrinsic orbital magnetic moment in a state with structure close to a uniform $`d_{x^2y^2}+id_{xy}`$. They furthermore showed, within a GL framework with coefficients determined by BCS weak-coupling theory, that this term is proportional to the particle-hole asymmetry of the normal metal from which the superconductor condenses, and dominates sufficiently far from the vortex core.

Interest in the possibility of order-parameter mixing in the vortex state was heightened by the experimental observation  of a plateau in the thermal conductivity as a function of the magnetic field $`H`$ when $`H`$ is above some critical value $`H^{}`$. Krishana et al., in particular, speculated that their observation of a sharp kink at $`H^{}`$ might be explained by the sudden onset of an out-of-phase $`d_{xy}`$ component at this critical field. The new high-field state was proposed to be fully gapped, with vanishing quasiparticle transport. There are several difficulties with this explanation, which we discuss below, but theorists were nonetheless persuaded to revisit the problem.

Laughlin  then pointed out, in analogy to the quantum Hall state, the peculiar nature of the time-reversal symmetry breaking $`d_{x^2y^2}+id_{xy}`$ state, which appears to be quite different from possible ground states in the case of the $`d,s`$-wave mixture. He proposed that the development of a magnetic moment coupling to the magnetic field might account for the phase transition, and put forward a $`T`$\- and $`H`$-dependent free energy functional driving the phase transition, and found a critical field $`H^{}T^2`$ similar to experiment. This special free energy functional does not have a microscopic basis and it is not clear yet whether there is such a field-induced secondary phase transition. To date, no magnetic-field induced transition has been found in relevant numerical studies of the vortex lattice.

Laughlin’s argument ignored the physics of the core region, thought to be negligible, but Ramakrishnan pointed out that quite similar effects are to be expected due to a combination of superfluid velocity Doppler shifts and Andreev reflection near the vortex cores. In higher fields, he proposed that these local $`d_{x^2y^2}+id_{xy}`$ patches might overlap, causing a transition to a uniform gapped state. This scenario is similar to one proposed by Movshovich et al. in the related case of magnetic impurities in a $`d`$-wave superconductor in zero field. Finally, Balatsky  has recently investigated the effect of the orbital Zeeman term near the upper critical field, and argued that the $`d`$-wave state is always unstable to a $`d_{x^2y^2}+id_{xy}`$ mixture. An unusual collective clapping mode in association with the relative phase of $`d_{xy}`$ to $`d_{x^2y^2}`$ was predicted in this superconductor.

Several important questions have not been addressed in the analyses of these issues thus far and have motivated this work. First and foremost, we would like to understand whether a phase transition of the type proposed by Krishana et al. is possible. Analyses in the GL regime are not applicable, and numerical calculations are not always useful to understand competing physical effects. We have therefore developed a systematic calculational approach capable of treating carefully the relevant quasiparticle states in the presence of spatially varying superflow together with the relevant subdominant order parameter components on an equal footing in the low-temperature phase. Our theory works for $`HH_{c2}`$, in which case the vortex core region can be safely neglected. Secondly, we would like to understand the structure of the vortex state and the role of the quasiparticles as a preliminary to the yet-unsolved problem of quasiparticle transport in applied magnetic field. In addition to being inapplicable at low temperatures, the GL calculations on which most of one’s intuition for this problem is based are unable to predict magnitudes of physical effects since they are based entirely on symmetry considerations. This is particularly important in the case of the orbital Zeeman coupling in the $`d_{x^2y^2},d_{xy}`$ mixing problem. The magnitude of the induced orbital moment is a very difficult quantity to estimate properly, as one might a priori deduce by analogy to the intrinsic orbital angular momentum problem in the $`{}_{}{}^{3}He`$ A-phase. In this case, it was found that naive calculations dramatically overestimated this effect, and we show here in fact that the orbital Zeeman coupling in the current problem is quite small.

In this paper, we adopt a semiclassical approach, expanding the BCS free energy in powers of the local superfluid velocity, local order parameter magnitude fluctuations, and their gradients. We thus neglect states possibly localized in the vortex core, and other quasiparticle bandstructure effects in a periodic vortex lattice discussed recently by several authors . Initially Volovik  and Kopnin and Volovik  proposed that a semiclassical analysis of this type should be valid only down to a scale $`(\mathrm{\Delta }_0^2/E_F)\sqrt{H/H_{c2}}`$. Recent numerical work  indicated, however, that the true crossover scale is much smaller for realistic systems with $`\mathrm{\Delta }_0/E_F1`$. Such fine details of the true quantum quasiparticle band structure will also be smeared out by impurity effects. We believe, therefore, that our neglect of the vortex core and quasiparticle bandstructure will be justified for cuprate superconductors at low fields ($`HH_{c2}`$) and temperatures, and that the current analysis will thus be adequate.

We begin by presenting in Sec. II the method, which involves a functional integral representation of the BCS free energy $`F`$, which we then expand in powers of slow superfluid velocity gradients and small subdominant order parameter components. Analytical results for $`F`$ in the GL regime, the low temperature regime, and, for $`s`$-wave case, an ultralow temperature regime where nonlinear superflow effects dominate, are given. This allows us in Sec. III to calculate the order parameter fluctuations directly. We then apply these results to the comparison of structure of a single isolated vortex with $`s`$ or $`d_{xy}`$ subdominant pairing at various temperatures in Sec. IV, and go on in Sec. V to discuss the prospects for observing a low temperature field-induced transition of this structure. In Sec. VI we discuss existing experiments and make some comments on the various available scenarios. In Appendix A, a detailed derivation of the free energy is presented, while Appendix B is devoted to a general calculation of the spontaneous magnetization in a $`d_{x^2y^2}+id_{xy}`$-wave superconductor.

## II FREE ENERGY

We start from a two-dimensional (2D) phenomenological BCS mean-field Hamiltonian in the mixed state:

$`\widehat{H}_{\mathrm{MF}}`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}{\displaystyle }d^2𝐫c_\sigma ^{}(𝐫)\{{\displaystyle \frac{1}{2m}}[i\stackrel{}{}{\displaystyle \frac{e}{c}}𝐀(𝐫)]^2\mu \}c_\sigma (𝐫)+{\displaystyle }d^2𝐫d^2𝐫^{}[\mathrm{\Delta }(𝐫,𝐫^{})c_{}^{}(𝐫)c_{}^{}(𝐫^{})+h.c.]`$ (2)
$`{\displaystyle d^2𝐫d^2𝐫^{}V(𝐫𝐫^{})|b(𝐫,𝐫^{})|^2},`$
where

$`V(𝐫)=V_d\mathrm{\Phi }_d(𝐫)+V_s+V_{xy}\mathrm{\Phi }_{xy}(𝐫)`$ (3)
with $`\mathrm{\Phi }_i`$ characterizing the irreducible representations $`d_{x^2y^2}`$, $`s`$, and $`d_{xy}`$, for which the lowest order basis functions over a circular Fermi surface are

$`\mathrm{\Phi }_{i𝐤}=\{\begin{array}{cc}\mathrm{cos}2\phi ,\hfill & i=d_{x^2y^2}\hfill \\ 1,\hfill & i=s\hfill \\ \mathrm{sin}2\phi ,\hfill & i=d_{xy}\hfill \end{array},`$ (7)
$`\mathrm{\Delta }(𝐫,𝐫^{})=V(𝐫𝐫^{})b(𝐫,𝐫^{})`$ is the pairing order parameter with $`b(𝐫,𝐫^{})=c_{}(𝐫^{})c_{}(𝐫)`$, and $`𝐀(𝐫)`$ the magnetic vector potential. Throughout the paper $`\mathrm{}=k_B=1`$ units are chosen. The magnetic field in the problem is perpendicular to the 2D plane, i.e., along the $`\widehat{z}`$ direction. We assume that $`V_d<0`$, and $`V_d,V_s`$ and $`V_{xy}`$ take such values that in the absence of the magnetic field, the superconducting state is of $`d_{x^2y^2}`$-wave symmetry with $`\mathrm{\Delta }(𝐫,𝐫^{})=\mathrm{\Delta }_0\mathrm{\Phi }_d(𝐫𝐫^{})`$. In the mixed state, $`\mathrm{\Delta }(𝐫,𝐫^{})`$ takes the form of $`\mathrm{\Delta }(𝐫,𝐫^{})=e^{i\varphi (\frac{𝐫+𝐫^{}}{2})}\stackrel{~}{\mathrm{\Delta }}(\frac{𝐫+𝐫^{}}{2},𝐫𝐫^{})`$, where $`\stackrel{~}{\mathrm{\Delta }}(𝐑,\stackrel{}{\rho })=\overline{\mathrm{\Delta }}_d(𝐑)\mathrm{\Phi }_d(\stackrel{}{\rho })+𝒟_s(𝐑)+𝒟_{xy}(𝐑)\mathrm{\Phi }_{xy}(\stackrel{}{\rho })`$, with $`\overline{\mathrm{\Delta }}_d(𝐑)=\mathrm{\Delta }_0+𝒟_d(𝐑)`$. The $`𝒟_i(𝐫)`$ are the magnetic field-induced pairing order parameter deviations from their values in zero field.

The partition function of Hamiltonian (2), after making a canonical transformation  to eliminate the phase field of the pairing order parameter $`\varphi (𝐫)`$,  is:

$`𝒵=𝒵_0\mathrm{exp}\left\{{\displaystyle \underset{n}{}}\mathrm{Tr}\mathrm{ln}\left[\widehat{M}_0+\left(\begin{array}{cc}\widehat{V}_1& \widehat{𝒟}\\ \widehat{𝒟}^{}& \widehat{V}_2\end{array}\right)\right]\right\},`$ (10)
where

$`𝒵_0`$ $`=`$ $`\mathrm{exp}(F_0/T),`$ (11)
$`F_0`$ $`=`$ $`{\displaystyle d^2𝐫\left(\frac{|\overline{\mathrm{\Delta }}_d(𝐫)|^2}{V_d}+\frac{|𝒟_s(𝐫)|^2}{V_s}+\frac{|𝒟_{xy}(𝐫)|^2}{V_{xy}}\right)},`$ (12)
and, in the momentum representation,

$`(\widehat{M}_0)_{𝐤,𝐤^{}}=\left(\begin{array}{cc}i\omega _n+ϵ_𝐤& \mathrm{\Delta }_𝐤\\ \mathrm{\Delta }_𝐤& i\omega _nϵ_𝐤\end{array}\right)\delta _{𝐤,𝐤^{}}`$ (15)
$`(\widehat{V}_1)_{𝐤,𝐤^{}}`$ $`=`$ $`{\displaystyle d^2𝐫e^{i(𝐤^{}𝐤)𝐫}\left[𝐯_s(𝐫)𝐤^{}+\frac{1}{2}mv_s^2(𝐫)\right]},`$ (16)
$`(\widehat{V}_2)_{𝐤,𝐤^{}}`$ $`=`$ $`{\displaystyle d^2𝐫e^{i(𝐤^{}𝐤)𝐫}\left[𝐯_s(𝐫)𝐤\frac{1}{2}mv_s^2(𝐫)\right]},`$ (17)
$`\widehat{𝒟}_{𝐤,𝐤^{}}`$ $`=`$ $`{\displaystyle \underset{i}{}}(\widehat{𝒟}_i)_{𝐤,𝐤^{}}={\displaystyle \underset{i}{}}{\displaystyle d^2𝐫e^{i(𝐤^{}𝐤)𝐫}𝒟_i(𝐫)\mathrm{\Phi }_{i\frac{𝐤+𝐤^{}}{2}}}`$ (18)
with $`\omega _n=(2n+1)\pi T`$ the fermion Matsubara frequency, $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }_0\mathrm{\Phi }_{d𝐤}`$, and

$`𝐯_s(𝐫)={\displaystyle \frac{e𝐀(𝐫)}{mc}}{\displaystyle \frac{\varphi (𝐫)}{2m}}`$ (19)
the supercurrent velocity. In writing down Eq. (10) we have used the relation $`\stackrel{}{}𝐯_s(𝐫)=0`$ which corresponds to the conservation of the supercurrent. Otherwise, $`mv_s^2(𝐫)/2`$ in Eqs. (16) and (17) has no significant effect and will be neglected hereafter.

The free energy resulting from Eq. (10) is

$`F`$ $`=`$ $`T\mathrm{Tr}\mathrm{ln}𝒵=F_0T\mathrm{Tr}\mathrm{ln}\widehat{M}_0+T{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{m}}\mathrm{Tr}\widehat{p}^m,`$ (20)
where

$`\widehat{p}=\widehat{g}\left(\begin{array}{cc}\widehat{V}_1& \widehat{𝒟}\\ \widehat{𝒟}^{}& \widehat{V}_2\end{array}\right),`$ (23)
with $`\widehat{g}`$ the Green-function matrix which, in the momentum representation, is

$`\widehat{g}_𝐤=\left(\begin{array}{cc}g_{1𝐤}& g_{2𝐤}\\ g_{2𝐤}& g_{4𝐤}\end{array}\right)=\left(\begin{array}{cc}\frac{i\omega _n+ϵ_𝐤}{W_{n𝐤}}& \frac{\mathrm{\Delta }_𝐤}{W_{n𝐤}}\\ \frac{\mathrm{\Delta }_𝐤}{W_{n𝐤}}& \frac{i\omega _nϵ_𝐤}{W_{n𝐤}}\end{array}\right).`$ (28)
Here,

$`W_{n𝐤}=\omega _n^2+E_𝐤^2,E_𝐤=\sqrt{ϵ_𝐤^2+\mathrm{\Delta }_k^2}.`$ (29)
The calculation of the trace of $`\widehat{p}^m`$ in Eq. (20) can be done by noting that

$`\mathrm{Tr}\widehat{p}^m`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{m}{}}}{\displaystyle \frac{d^2𝐤_j}{(2\pi )^2}d^2𝐫_je^{i(𝐤_j𝐤_{j+1})𝐫_j}\mathrm{Tr}\left(\underset{j=1}{\overset{m}{}}\widehat{p}_{𝐤_j,𝐤_{j+1}}(𝐫_j)\right)},`$ (30)
where $`𝐤_{m+1}=𝐤_1`$. In most of the bulk region not close to a vortex core, $`𝐯_s`$ as well as the $`𝐯_s`$-induced $`𝒟_i`$ are spatially slowly varying functions. Thus we are allowed to expand $`𝐯_s(𝐫_j)`$ and $`𝒟_i(𝐫_j)`$ in $`\widehat{p}_{𝐤,𝐤_1}(𝐫)`$ in Eq. (30) as power series in their derivatives ,

$`𝐯_s(𝐫_j)`$ $``$ $`e^{(𝐫_j𝐫)\stackrel{}{}_𝐫}𝐯_s(𝐫)=𝐯_s(𝐫)+[(𝐫_j𝐫)\stackrel{}{}_𝐫]𝐯_s(𝐫)+\mathrm{},`$ (31)
$`𝒟_i(𝐫_j)`$ $``$ $`e^{(𝐫_j𝐫)\stackrel{}{}_𝐫}𝒟_i(𝐫)=𝒟_i(𝐫)+[(𝐫_j𝐫)\stackrel{}{}_𝐫]𝒟_i(𝐫)+\mathrm{},`$ (32)
the first few terms of which make main contribution to Tr$`\widehat{p}^m`$. This property enables us to develop a perturbation theory to obtain the free energy. The resulting calculation is straightforward but tedious, and is summarized in Appendix A. The final result for the free energy with respect to the pairing order parameters $`\stackrel{~}{F}=d^2𝐫\stackrel{~}{f}(𝐫)`$, where the free energy density $`\stackrel{~}{f}(𝐫)`$, keeping terms up to quadratic in $`𝒟_i`$, is

$`\stackrel{~}{f}(𝐫)=f_{\overline{\mathrm{\Delta }}_d}(𝐫)+f_s(𝐫)+f_{xy}(𝐫)+\delta f(𝐫),`$ (33)
$`f_{\overline{\mathrm{\Delta }}_d}(𝐫)={\displaystyle \frac{|\overline{\mathrm{\Delta }}_d(𝐫)|^2}{V_d}}T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\mathrm{ln}\left(\frac{\stackrel{~}{W}_{n𝐤}(𝐫)}{W_{n𝐤}}\right)},`$ (34)
$`f_s(𝐫)=\mathrm{\Delta }_0^2(T)N_0\{L_s^{}(T,𝐯_s)\overline{𝒟}_s^{}(𝐫)+L_s^{\prime \prime }(T,𝐯_s)\overline{𝒟}_s^{\prime \prime }(𝐫)+[c_s+\eta _s^{}(T,𝐯_s)][\overline{𝒟}_s^{}(𝐫)]^2+[c_s+\eta _s^{\prime \prime }(T,𝐯_s)][\overline{𝒟}_s^{\prime \prime }(𝐫)]^2\},`$ (35)
$`f_{xy}(𝐫)=F_4(𝐫)+F^{OZ}(𝐫)+\mathrm{\Delta }_0^2(T)N_0\left\{[c_{xy}+\eta _{xy}^{}(T,𝐯_s)][\overline{𝒟}_{xy}^{}(𝐫)]^2+[c_{xy}+\eta _{xy}^{\prime \prime }(T,𝐯_s)][\overline{𝒟}_{xy}^{\prime \prime }(𝐫)]^2\right\},`$ (36)
$`F^{OZ}(𝐫)=\mathrm{\Delta }_0^2(T)N_0Q^{OZ}(T){\displaystyle \frac{e}{mc}}B(𝐫)\overline{𝒟}_{xy}^{\prime \prime }(𝐫),`$ (37)
$`F_4(𝐫)=\mathrm{\Delta }_0^2(T)N_0\left\{L_{xy}^{}(T,𝐯_s)\overline{𝒟}_{xy}^{}(𝐫)+L_{xy}^{\prime \prime }(T,𝐯_s)\overline{𝒟}_{xy}^{\prime \prime }(𝐫)\right\},`$ (38)
where  and <sup>′′</sup> indicate real and imaginary parts, respectively, $`\overline{𝒟}_i=𝒟_i/\mathrm{\Delta }_0(T)`$, $`i=s,d_{xy}`$, $`c_s=(V_sN_0)^12c_d`$, $`c_{xy}=(V_{xy}N_0)^1c_d`$, with $`c_d=(V_dN_0)^1`$, $`B(𝐫)`$ is the magnetic induction, and $`\stackrel{~}{W}_{n𝐤}(𝐫)`$, $`Q^{OZ}(T)`$, $`L_i(T,𝐯_s)`$ and $`\eta _i(T,𝐯_s)`$ are defined in Appendix A. $`f_{\overline{\mathrm{\Delta }}_d}`$ in Eq. (34) is the free energy density in association with the dominant $`d_{x^2y^2}`$ component in the absence of $`𝒟_s`$ and $`𝒟_{xy}`$. $`f_s`$ and $`f_{xy}`$ are the free energy densities for the $`s`$ and $`d_{xy}`$ components, respectively. $`\delta f`$ involves terms of high orders in $`v_s`$, $`𝒟_i`$, derivatives of $`v_s`$, and mixed $`s`$ and $`d_{xy}`$ terms. Note in Eqs. (35) and (36) the quadratic terms in $`D_i`$, $`i=s,d_{xy}`$ can be easily reformulated to coincide with the familiar GL form: $`(c_i+\frac{\eta ^{}+\eta ^{\prime \prime }}{2})|\overline{𝒟}_i(𝐫)|^2+\frac{\eta ^{}\eta ^{\prime \prime }}{4}[\overline{𝒟}_i(𝐫)^2+\overline{𝒟}_i^{}(𝐫)^2]`$.  In addition, derivatives of $`𝒟_i`$ terms are absorbed into the powers of $`𝒟_i`$ terms by partial integration, as indicated in Appendix A.

With the free energy in Eqs. (33)-(38) we are in a position to investigate the vortex state. We will first show general results for an arbitrary $`𝐯_s(𝐫)`$ distribution  in Sec. III, and apply the theory to the single vortex case in Sec. IV.

## III Order parameters

### A $`d`$-wave pairing order parameter $`\overline{\mathrm{\Delta }}_d(𝐫)=\mathrm{\Delta }_0+𝒟_d(𝐫)`$

In studying the dominant $`d_{x^2y^2}`$ component, $`𝒟_s`$ and $`𝒟_{xy}`$ can be set to zero since the $`\overline{\mathrm{\Delta }}_d`$-$`𝒟_s`$ and $`\overline{\mathrm{\Delta }}_d`$-$`𝒟_{xy}`$ mixing terms appear in higher orders of $`v_s`$ or its derivatives. The gap equation $`f_{\overline{\mathrm{\Delta }}_d}/\overline{\mathrm{\Delta }}_d(𝐫)=0`$ produces

$`{\displaystyle \frac{\overline{\mathrm{\Delta }}_d(𝐫)}{V_d}}=T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\frac{\overline{\mathrm{\Delta }}_d(𝐫)\mathrm{cos}^22\phi }{\stackrel{~}{W}_{n𝐤}(𝐫)}},`$ (39)
where $`\stackrel{~}{W}_{n𝐤}(𝐫)`$ is defined in Eq. (A6) in Appendix A. At $`H=0`$, Eq. (39) reduces to

$`1=V_dT{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\frac{\mathrm{cos}^22\phi }{W_{n𝐤}}}=V_d{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\frac{\mathrm{cos}^22\phi \mathrm{tanh}(E_𝐤/2T)}{2E_𝐤}},`$ (40)
which yields the well-known $`T_c`$ formula $`c_d=(1/2)\mathrm{ln}(2e^\gamma \omega _D/\pi T_c)`$, with $`\omega _D`$ the BCS cutoff and $`\gamma 0.5772`$ the Euler constant, as well as the asymptotic behaviors of the gap maximum

$$\mathrm{\Delta }_0(T)\{\begin{array}{cc}3.54[T_c(TT_c)]^{1/2},\hfill & T\mathrm{\Delta }_0(T)\hfill \\ 2.14T_c0.39(T/T_c)^3T_c,\hfill & T\mathrm{\Delta }_0(T)\hfill \end{array}.$$
(41)
The presence of vortices depletes the pairing order parameter $`\mathrm{\Delta }_0`$ by $`𝒟_d(𝐫)`$. Sufficiently far from the vortex cores, either $`T`$ or $`\mathrm{\Delta }_0(T)`$ is larger than $`v_s(𝐫)k_F`$ in the whole $`T`$ regime, and thus we can treat $`v_s(𝐫)k_F`$ and $`𝒟_d(𝐫)`$ perturbatively. It is straightforward to show that the normalized real part of $`𝒟_d(𝐫)`$ is

$`{\displaystyle \frac{𝒟_d^{}(𝐫)}{\mathrm{\Delta }_0(T)}}\{\begin{array}{cc}\frac{1}{3/4+ϵ_Fmv_s^2(𝐫)/\mathrm{\Delta }_0^2(T)}\frac{ϵ_Fmv_s^2(𝐫)}{\mathrm{\Delta }_0^2(T)},\hfill & \mathrm{GL}\mathrm{regime}\hfill \\ \frac{2(\mathrm{ln}\mathrm{\hspace{0.17em}2})}{1[9\zeta (3)T^34(\mathrm{ln}\mathrm{\hspace{0.17em}2})ϵ_FTmv_s^2(𝐫)]/\mathrm{\Delta }_0^3(T)}\frac{ϵ_Fmv_s^2(𝐫)}{\mathrm{\Delta }_0^2(T)}\frac{T}{\mathrm{\Delta }_0(T)},\hfill & v_sk_FT\mathrm{\Delta }_0(T)\hfill \\ \frac{_{l=\pm 1}|\mathrm{cos}(\theta +l\frac{\pi }{4})|^3}{3[_{l=\pm 1}|\mathrm{cos}(\theta +l\frac{\pi }{4})|^3][v_s(𝐫)k_F]^3/\mathrm{\Delta }_0^3(T)}\frac{ϵ_Fmv_s^2(𝐫)}{\mathrm{\Delta }_0^2(T)}\frac{v_sk_F}{\mathrm{\Delta }_0(T)},\hfill & Tv_sk_F\mathrm{\Delta }_0(T)\hfill \end{array},`$ (45)
where $`ϵ_F`$ is the Fermi energy and $`\zeta (3)`$ is the Riemann function. The imaginary part of $`𝒟_d(𝐫)`$ driven by the derivatives of the supercurrent is unimportant compared with $`𝒟_d^{}(𝐫)`$ in the whole temperature region, for its driving term $`F_d^{(1)}`$ in Eq. (A16) is very small as discussed in Appendix A.

### B Field-Induced $`𝒟_s(𝐫)`$

From Eqs. (33) and (35) we obtain the gap equation for the $`s`$ component, $`f_s(𝐫)/𝒟_s(𝐫)=0`$, which gives

$`\overline{𝒟}_s^{}(𝐫){\displaystyle \frac{L_s^{}(T,𝐯_s)}{2[c_s+\eta _s^{}(T,𝐯_s)]}},\overline{𝒟}_s^{\prime \prime }(𝐫){\displaystyle \frac{L_s^{\prime \prime }(T,𝐯_s)}{2[c_s+\eta _s^{\prime \prime }(T,𝐯_s)]}},`$ (46)
where $`L_s^{}(T,𝐯_s)`$, $`L_s^{\prime \prime }(T,𝐯_s)`$, $`\eta _s^{}(T,𝐯_s)`$ and $`\eta _s^{\prime \prime }(T,𝐯_s)`$ are defined in Appendix A.

For $`T>𝐯_s𝐤_F`$, we can plug $`L_s`$ and $`\eta _s`$ obtained in Eqs. (A18)-(A21) into Eq. (46) to get the scaling functions of $`\overline{𝒟}_s^{}(𝐫)`$ and $`\overline{𝒟}_s^{\prime \prime }(𝐫)`$,

$`\overline{𝒟}_s^{}(𝐫)`$ $``$ $`G_{s1}({\displaystyle \frac{\mathrm{\Delta }_0(T)}{T}},{\displaystyle \frac{ϵ_Fmv_s^2(𝐫)}{T^2}}){\displaystyle \frac{ϵ_Fm[v_{sx}^2(𝐫)v_{sy}^2(𝐫)]}{T^2}},`$ (47)
$`\overline{𝒟}_s^{\prime \prime }(𝐫)`$ $``$ $`G_{s2}({\displaystyle \frac{\mathrm{\Delta }_0(T)}{T}},{\displaystyle \frac{ϵ_Fmv_s^2(𝐫)}{T^2}}){\displaystyle \frac{ϵ_F_xv_{sx}(𝐫)}{T^2}},`$ (48)
where

$`G_{s1}(d,z)`$ $`=`$ $`{\displaystyle \frac{h_{1s}(d)}{c_sh_{3s}(d)+2d^2h_{2s}(d)+zh_{4s}(d)}},`$ (49)
$`G_{s2}(d,z)`$ $`=`$ $`{\displaystyle \frac{2h_{2s}(d)}{c_sh_{3s}(d)+zh_{4s}(d)}},`$ (50)
with $`h_{1s}`$, $`\mathrm{}`$, $`h_{4s}`$ defined in Eqs. (A22) and (A23).

It is easy to see that in the GL regime,

$`\overline{𝒟}_s^{\mathrm{GL}}(𝐫){\displaystyle \frac{1}{c_s+\frac{3}{4}\frac{\mathrm{\Delta }_0^2(T)}{\pi ^2T^2}+2\frac{ϵ_Fmv_s^2(𝐫)}{\pi ^2T^2}}}{\displaystyle \frac{ϵ_Fm[v_{sx}^2(𝐫)v_{sy}^2(𝐫)]}{\pi ^2T^2}},`$ (51)
$`\overline{𝒟}_s^{\prime \prime \mathrm{GL}}(𝐫){\displaystyle \frac{1}{c_s\frac{\mathrm{\Delta }_0^2(T)}{4\pi ^2T^2}+2\frac{ϵ_Fmv_s^2(𝐫)}{\pi ^2T^2}}}{\displaystyle \frac{ϵ_F_xv_{sx}(𝐫)}{\pi ^2T^2}}.`$ (52)
It follows that the ratio between the real and imaginary parts of $`𝒟_s`$

$`{\displaystyle \frac{\overline{𝒟}_s^{\mathrm{GL}}(𝐫)}{\overline{𝒟}_s^{\prime \prime \mathrm{GL}}(𝐫)}}{\displaystyle \frac{m[v_{sx}^2(𝐫)v_{sy}^2(𝐫)]}{_xv_{sx}(𝐫)}}`$ (53)
is of order unity in most of the region outside the core where $`v_s(𝐫)1/r_0`$ and $`_xv_{sx}(𝐫)1/r_0^2`$ with $`r_0`$ the distance to the closest vortex core center.

In the case of $`𝐯_s𝐤_F<T\mathrm{\Delta }_0(T)`$,

$`\overline{𝒟}_s^{\mathrm{lowT}}(𝐫){\displaystyle \frac{(\mathrm{ln}\mathrm{\hspace{0.17em}2})}{c_s+\frac{1}{2}+\frac{ϵ_Fmv_s^2(𝐫)}{4T\mathrm{\Delta }_0}}}{\displaystyle \frac{T}{\mathrm{\Delta }_0}}{\displaystyle \frac{ϵ_Fm[v_{sx}^2(𝐫)v_{sy}^2(𝐫)]}{\mathrm{\Delta }_0^2}},`$ (54)
$`\overline{𝒟}_s^{\prime \prime \mathrm{lowT}}(𝐫){\displaystyle \frac{12(\mathrm{ln}\mathrm{\hspace{0.17em}2})\frac{T}{\mathrm{\Delta }_0}}{c_s\frac{1}{2}+2(\mathrm{ln}\mathrm{\hspace{0.17em}2})\frac{T}{\mathrm{\Delta }_0}+\frac{ϵ_Fmv_s^2(T)}{4T\mathrm{\Delta }_0}}}{\displaystyle \frac{ϵ_F_xv_{sx}(𝐫)}{\mathrm{\Delta }_0^2}}.`$ (55)
So the ratio between the real and imaginary parts becomes

$`{\displaystyle \frac{\overline{𝒟}_s^{\mathrm{lowT}}(𝐫)}{\overline{𝒟}_s^{\prime \prime \mathrm{lowT}}(𝐫)}}(\mathrm{ln}\mathrm{\hspace{0.17em}2}){\displaystyle \frac{T}{\mathrm{\Delta }_0}}{\displaystyle \frac{m[v_{sx}^2(𝐫)v_{sy}^2(𝐫)]}{_xv_{sx}(𝐫)}}.`$ (56)
Comparing Eqs. (53) and (56) we find an extra prefactor $`T/\mathrm{\Delta }_0`$ is acquired at low $`T`$, indicative of a suppressed real part of the $`s`$ component with decreasing $`T`$. This is an interesting observation in the present work, the consequence of which on the structure of a single vortex will be discussed in Sec. IV.

For $`T<𝐯_s𝐤_F`$, which can be achieved either by lowering $`T`$ in a certain spatial position or approaching the core region at a certain $`T`$, the nonlinear effects dominate over the thermal effects and the prefactor $`T/\mathrm{\Delta }_0`$ in Eq. (56) is expected to be replaced by $`v_sk_F/\mathrm{\Delta }_0`$. So the real part is negligibly small, and the imaginary part at $`T=0`$ is

$`\overline{𝒟}_s^{\prime \prime T=0}(𝐫){\displaystyle \frac{1S_\theta v_s(𝐫)k_F/\mathrm{\Delta }_0}{c_s1/2+S_\theta v_s(𝐫)k_F/2\mathrm{\Delta }_0}}{\displaystyle \frac{ϵ_F_xv_{sx}(𝐫)}{\mathrm{\Delta }_0^2}},`$ (57)
where $`S_\theta =_{l=\pm 1}\left|\mathrm{cos}(\theta +l\frac{\pi }{4})\right|`$. It is interesting to find from Eqs. (57) and (55) that the factor $`T/\mathrm{\Delta }_0`$ at $`𝐯_s𝐤_F<T`$ is replaced by $`v_sk_F/\mathrm{\Delta }_0`$ (up to some prefactor) at $`T<𝐯_s𝐤_F`$, reflecting the nonlinear effect due to the Doppler energy shift. Its effect on the free energy and the penetration depth in the Meissner state of a $`d`$-wave superconductor was extensively discussed in the previous papers by the present authors.

### C Field-Induced $`𝒟_{xy}(𝐫)`$

As shown in Eq. (36), there are two terms, $`F^{OZ}`$ and $`F_4`$, competing in driving the induced $`d_{xy}`$ component. The first term $`F^{OZ}`$, the so-called orbital Zeeman term, can be rewritten as

$$F^{OZ}(𝐫)=𝐌(𝐫)𝐁(𝐫),$$
(58)
with $`𝐌(𝐫)`$ the effective magnetic moment associated with the internal orbital current of Cooper pairs in a $`d_{x^2y^2}+id_{xy}`$-wave superconductor. A general derivation of this spontaneous magnetization is presented in Appendix B. As also discussed there, $`𝐌(𝐫)`$ is proportional to the particle-hole asymmetry $`\alpha `$ (see Eq. (A44)), making the Zeeman orbital effects very small for general density of states. It is therefore worth pointing out that analyses of the $`d_{x^2y^2}`$-$`d_{xy}`$ mixing problem based solely on a symmetry-based analysis of the orbital moment  may lead to unrealistic results. Besides, it is important to note the very weak temperature dependence of the coefficient $`Q^{OZ}`$ of $`F^{OZ}`$, as shown in Fig. 1.

The second term, $`F_4`$, is in the order $`𝒪(v_s^4)`$ and contains driving terms for $`𝒟_{xy}^{\prime \prime }`$ as well as $`𝒟_{xy}^{}`$. This term can be significant particularly if the system has small or zero particle-hole asymmetry. Note that in the $`s`$-wave case, it is irrelevant except for very short length scales of order the core size, due to the nonvanishing leading ($`𝒪(v_s^2)`$) term.

The orbital Zeeman term has been invoked by Laughlin  and Balatsky  as the effect driving a putative transition to a time reversal symmetry-breaking $`d_{x^2y^2}+id_{xy}`$ state induced by field. It is therefore particularly interesting to investigate the relative importance of $`F^{OZ}`$ and $`F_4`$ within the current approach. For fields $`HH_{c1}`$, the overlap of vortices leads to nearly homogeneous $`B(𝐫)H`$ in space and hence $`F^{OZ}/|𝒟_{xy}|`$ can be taken as a constant over the bulk. On the other hand, $`F_4/|𝒟_{xy}|`$, which scales as $`(\xi /r)^4`$, is strongly space dependent. It increases rapidly when approaching the vortex core, but decays into the bulk. Thus we expect a critical radius $`r^{}`$ beyond which $`F^{OZ}`$ dominates over $`F_4`$, but for $`r<r^{}`$ $`F_4`$ becomes more important and determines the structure. $`r^{}`$ can be estimated from $`F_4(r^{})/F^{OZ}(r^{})1`$. It turns out to grow with increasing $`H`$ and decreasing $`T`$. In the GL and low $`T`$ regime, we have

$$r^{}(T,H)=\{\begin{array}{cc}\xi \left(\frac{ϵ_F}{eH/mc}\right)^{1/4}\xi \left(\frac{R_H}{\lambda _F}\right)^{1/2},\hfill & GL\hfill \\ \xi \left(\frac{\mathrm{\Delta }_0(T)}{T}\right)^{1/4}\left(\frac{ϵ_F}{eH/mc}\right)^{1/4}\xi \left(\frac{\mathrm{\Delta }_0(T)}{T}\right)^{1/4}\left(\frac{R_H}{\lambda _F}\right)^{1/2},\hfill & \mathrm{low}T\hfill \end{array},$$
(59)
with $`\lambda _F`$ the Fermi wavelength, $`\xi =v_F/\pi \mathrm{\Delta }_0`$ the superconducting coherence length, and $`R_H=\sqrt{c/eH}`$ the average intervortex distance. We summarize the above estimates in Table I, and display schematically the competition between the two driving terms in Fig. 2. Note the “low T” results for $`F_4`$ obtain only down to temperatures above the local superfluid velocity Doppler shift $`v_sk_F`$. For lower $`T`$, the perturbation calculation for $`F_4`$ that we present in Appendix A.3 breaks down because of the zero modes in $`E_𝐤`$. Instead one has to get a full expression for the local Doppler shift and the derivatives of the supercurrents which is suitable for perturbation expansion after integrating over momenta. This is apparently too complicated to be achieved analytically and thus requires a self-consistent numerical work, which is beyond the scope of the present paper. However, based on the nonlinear results for the $`s`$ component and $`Q^{OZ}`$, we do not expect anything qualitatively new in the extremely low $`T`$ case compared with the $`𝐯_s𝐤_FT\mathrm{\Delta }_0(T)`$ case.

The resulting $`d_{xy}`$-wave order parameter, for $`\xi <r<r^{}`$ where $`F_4`$ dominates, is

$`\overline{𝒟}_{xy}^{}(𝐫){\displaystyle \frac{L_{xy}^{}(T,𝐯_s)}{2[c_{xy}+\eta _{xy}^{}(T,𝐯_s)]}},\overline{𝒟}_{xy}^{\prime \prime }(𝐫){\displaystyle \frac{L_{xy}^{\prime \prime }(T,𝐯_s)}{2[c_{xy}+\eta _{xy}^{\prime \prime }(T,𝐯_s)]}},`$ (60)
where $`L_{xy}^{}(T,𝐯_s)`$, $`L_{xy}^{\prime \prime }(T,𝐯_s)`$, $`\eta _{xy}^{}(T,𝐯_s)`$, and $`\eta _{xy}^{\prime \prime }(T,𝐯_s)`$ are defined in Eqs. (A41), (A42), (A45) and (A46), and their asmptotic behaviors can be found in Table II. In the region $`r^{}<r<R_H`$ where $`F^{OZ}`$ dominates, the $`d_{xy}`$ component reads

$`\overline{𝒟}_{xy}^{}(𝐫)0,\overline{𝒟}_{xy}^{\prime \prime }(𝐫){\displaystyle \frac{Q^{OZ}(T)}{2[c_{xy}+\eta _{xy}^{\prime \prime }(T,𝐯_s)]}}{\displaystyle \frac{eH}{mc}},`$ (61)
where $`Q^{OZ}(T)`$ defined in Eq. (A43). Eq. (61) and the asymptotic behaviors of $`Q^{OZ}`$, $`\eta _{xy}^{}(T,𝐯_s)`$, and $`\eta _{xy}^{\prime \prime }(T,𝐯_s)`$ shown in Tabel II lead to

$`\overline{𝒟}_{xy}^{\prime \prime }(𝐫)\{\begin{array}{cc}\frac{1}{2}\alpha (2c_d1)[c_{xy}\frac{3}{8}\frac{\mathrm{\Delta }_0^2(T)}{\pi ^2T^2}+\frac{ϵ_Fmv_s^2(𝐫)}{\pi ^2T^2}]^1\frac{1}{ϵ_F}\frac{eH}{mc}\hfill & \mathrm{\Delta }_0(T)T\hfill \\ \frac{1}{2}\alpha [\mathrm{ln}(\frac{4\omega _D}{\mathrm{\Delta }_0})\frac{1}{2}6(\mathrm{ln2})\frac{T}{\mathrm{\Delta }_0}][c_{xy}\frac{1}{2}+2(\mathrm{ln2})\frac{T}{\mathrm{\Delta }_0}+\frac{ϵ_Fmv_s^2(𝐫)}{4\mathrm{\Delta }_0T}]^1\frac{1}{ϵ_F}\frac{eH}{mc}\hfill & T\mathrm{\Delta }_0(T)\hfill \end{array}.`$ (64)
It is interesting to estimate the value of $`r^{}`$ in a real material. It follows from Eq. (59) that for $`TT_c`$, $`r^{GL}/R_H\sqrt{\xi /\lambda _F}(H/H_{c2})^{1/4}`$. For high-$`T_c`$ cuprates, $`\xi /\lambda _F10`$, and thus $`r^{GL}>R_H`$ for fields $`H>0.01H_{c2}`$. In materials with larger $`\xi /\lambda _F`$, $`r^{GL}`$ becomes order of $`R_H`$ for smaller fields. Since from Eq. (59) $`r^{}`$ increases with decreasing temperatures, it seems unlikely that the orbital Zeeman free energy plays an important role in determining the local order parameter in the vortex state for fields in the Tesla range.

We would like to make some further remarks on the magnetic field dependence of $`\overline{𝒟}_{xy}^{\prime \prime }`$ for $`r>r^{}`$ at low temperatures. $`\overline{𝒟}_{xy}^{\prime \prime }`$ shown in Eqs. (61) and (64) are obtained from minimizing the free energy density $`f_{xy}(𝐫)`$ in Eq. (36). This free energy density is up to quadratic in $`\overline{𝒟}_{xy}^{\prime \prime }`$, which is sufficient for low field and for generic $`c_{xy}1`$. In the case of larger field and/or special case of $`c_{xy}`$ close to or smaller than 1, one has to include the free energy density term cubed in $`\overline{𝒟}_{xy}^{\prime \prime }`$. This term can be easily found to be $`(\overline{𝒟}_{xy}^{\prime \prime })^3/3`$ with the spatial dependence of the coefficient neglected. Thus the free energy density up to cubed in $`\overline{𝒟}_{xy}^{\prime \prime }`$ for $`r>r^{}`$ at low temperatures can be written as

$`\stackrel{~}{f}_{xy}=\gamma B\overline{𝒟}_{xy}^{\prime \prime }+c_0(\overline{𝒟}_{xy}^{\prime \prime })^2+{\displaystyle \frac{1}{3}}(\overline{𝒟}_{xy}^{\prime \prime })^3,`$ (65)
where $`c_0=c_{xy}+\eta _{xy}^{\prime \prime }(T,𝐯_s)`$ and $`\gamma =\mathrm{\Delta }_0^2N_0Q^{OZ}e/mc`$. Minimizing $`\stackrel{~}{f}_{xy}`$ with respect to $`\overline{𝒟}_{xy}^{\prime \prime }`$ we immediately get

$`\overline{𝒟}_{xy}^{\prime \prime }\sqrt{c_0^2+\gamma B/2}c_0.`$ (66)
It is obvious that there is crossover of the linear-$`B`$ dependence of $`\overline{𝒟}_{xy}^{\prime \prime }`$ for $`\gamma B/2c_0^2`$ to square root of $`B`$ dependence of $`\overline{𝒟}_{xy}^{\prime \prime }`$ for $`\gamma B/2c_0^2`$. This interesting behavior is shown in Fig. (3).

In a periodic vortex lattice, the supercurrent field $`𝐯_s(𝐫)`$ in the London approximation reads

$`𝐯_s(𝐫)={\displaystyle \frac{\pi }{m}}{\displaystyle \underset{𝐊G}{}}{\displaystyle \frac{i\stackrel{}{}_𝐫(e^{i𝐊𝐫})}{K^2+\lambda ^2}},`$ (67)
and thus there will be special symmetry points where $`𝐯_s=0`$. At these points, a careful examination shows that the coefficients $`L_{xy}^{}`$ and $`L_{xy}^{\prime \prime }`$ in $`F_4`$ (Eq. (38)) vanish, and that $`𝒟_{xy}`$ is driven entirely by $`F^{OZ}`$, as shown in Eq. (66). This coincides with the numerical result of Yasui and Kita.

Up to now, we have not included the effect of the $`d_{x^2y^2}`$ order-parameter suppression $`𝒟_d`$ on the $`d_{xy}`$ component, since we have shown in Sec. III A that $`𝒟_d`$ is negligibly small in the bulk region, in which we are primarily interested. However, the mechanism for a transition to a $`d_{x^2y^2}+id_{xy}`$ state proposed by Ramakrishnan involves precisely the interplay between $`𝒟_d`$ and the supercurrent near the core, leading potentially to local $`d_{x^2y^2}+id_{xy}`$ patches which then overlap at some critical field. Motivated by this suggestion, we examine the free energy terms including this effect within our approximation. These terms, which we refer to as the order parameter suppression terms $`F^{OPS}`$, are obtained from Eq. (A9-A12) (and have been already neglected in arriving at Eq. (20)). It can easily be seen that the $`𝒟_{xy}`$ component couples to derivatives of either $`𝒟_d(𝐫)`$ or of $`𝐯_s`$ in these terms, indicative of a pure nonlocal effect as found by Ramakrishnan. However, for $`𝐯_s𝐤_F<T`$, our surprising finding is that up to leading order, terms including derivatives of $`𝒟_d(𝐫)`$ vanish, leaving

$`F^{OPS}={\displaystyle d^2r\left(T\underset{n}{}\frac{d^2𝐤}{(2\pi )^2}\frac{ϵ_𝐤\mathrm{cos}^22\phi }{W_{n𝐤}^2}\right)𝒟_d(𝐫)[_𝐫\times 𝐯_s(𝐫)]_z𝒟_{xy}^{\prime \prime }(𝐫)}.`$ (68)
Since the $`[_𝐫\times 𝐯_s(𝐫)]_z`$ factor is dominated by the vector potential $`𝐀(𝐫)`$ rather than $`\varphi `$ part in $`𝐯_s(𝐫)`$, it may be replaced by $`eB(𝐫)/(mc)`$. Comparing Eqs. (68) and (37), we see that $`F^{OPS}`$ in the core region is of the same order as $`F^{OZ}`$ in the bulk (in particular, it is also proportional to the particle-hole asymmetry of the normal state), and may be viewed as the leading correction to $`F^{OZ}`$, if $`\mathrm{\Delta }_0(T)`$ in $`F^{OZ}`$ is replaced by $`\overline{\mathrm{\Delta }}_d(𝐫)=\mathrm{\Delta }_0(T)+𝒟_d(𝐫)`$. It is therefore clear from our previous discussion of $`F^{OZ}`$ that $`F^{OPS}`$ also gives in fact a very small effect even near the vortex core.

## IV Single vortex

The results obtained in Sec. III enable us to compare the subdominant order parameters, at various temperatures, in the presence of vortices characterized by a certain superfluid velocity field $`𝐯_s(𝐫)`$.  Our purpose in this section is to test these results in the concrete case of a single isolated vortex.

We use cylindrical coordinates $`𝐫=(r,\theta )`$ with the origin located at the vortex core center. The phase field of the order parameter $`\varphi (𝐫)=\theta `$ leading to $`\varphi (𝐫)=(1/2mr)\widehat{\theta }`$. In the spatial regime of interest, the magnetic field is roughly homogeneous and thus $`𝐀(𝐫)=Br\widehat{\theta }`$. For $`rR_H`$, $`|\varphi (𝐫)/2m|eA(𝐫)/(mc)`$, so we can neglect $`𝐀`$, and simply write the supercurrent components and their derivatives as

$$v_{sx}\frac{\mathrm{sin}\theta }{2mr},v_{sy}\frac{\mathrm{cos}\theta }{2mr},_xv_{sx}(𝐫)\frac{\mathrm{sin}2\theta }{2mr^2},$$
(69)
except for studying the orbital Zeeman term, in which case $`A(𝐫)`$ is no longer negligible because $`\times 𝐀`$ enters.

### A $`d_{x^2y^2}`$-$`s`$ mixing

In the GL regime, the $`s`$-wave subdominant order parameter has been investigated by many authors,  with substantial agreement. The existing results in this limit are easily shown to be recovered in the present theory. Inserting Eq. (69) into Eqs. (51) and (52) leads, in the generic case of $`c_s1`$, to

$$\overline{𝒟}_s^{GL}(𝐫)0.6c_s^1\left(\frac{\xi }{r}\right)^2(\mathrm{cos}2\theta +2i\mathrm{sin}2\theta )=0.3c_s^1\left(\frac{\xi }{r}\right)^2(3e^{2i\theta }e^{2i\theta })$$
(70)
which is consistent with the results in Ref. in the bulk asymptotic region $`r\xi `$. As a result, the $`s`$ component is of four-fold symmetry, and the relative winding of $`s`$ and $`d`$ components is uniform across the whole vortex as shown in Fig. 4.

In the low $`T`$ case, as discussed in Sec. III B, $`𝒟_s^{}(𝐫)`$ is smaller than $`𝒟_s^{\prime \prime }(𝐫)`$ by a factor of (ln 2)$`T/\mathrm{\Delta }_0`$, and the vortex structure at low $`T`$ is thus expected to be qualitatively different from that in the GL regime. From Eqs. (69), (54), and (55), one finds that

$$\overline{𝒟}_s^{\mathrm{lowT}}(𝐫)1.23c_s^1\left(\frac{\xi }{r}\right)^2\left[(\mathrm{ln2})\frac{T}{\mathrm{\Delta }_0}\mathrm{cos}2\theta +2i\mathrm{sin}2\theta \right].$$
(71)
The magnitude and phase of $`𝒟_s^{\mathrm{lowT}}`$ are shown in Fig. 4. The relative winding of $`s`$ and $`d`$ components at low $`T`$ takes place in a very narrow region of real space near antinode directions set by max$`(T/\mathrm{\Delta }_0,v_sk_F/\mathrm{\Delta }_0)`$. In Fig. 5, we show the $`T`$ dependence of $`\overline{𝒟}_s^{\prime \prime }`$ and $`\overline{𝒟}_s^{}`$ in a spatial position in the bulk.

### B $`d_{x^2y^2}`$-$`d_{xy}`$ mixing

As discussed in Sec. III C, the competition between $`F^{OZ}`$ and $`F_4`$ divides the outside-vortex-core region into two rings: the outer ring region $`r^{}<r<R_H`$ is dominated by $`F^{OZ}`$ resulting in a rigid $`d_{x^2y^2}+id_{xy}`$ superconducting state with spatially nearly constant $`𝒟_{xy}^{\prime \prime }`$ obtained in Eqs. (61) and (64); In the inner ring $`\xi <r<r^{}`$, $`F_4`$ is more important and a spatially varying $`𝒟_{xy}`$ is expected. Now we show the single vortex structure in the inner ring region. We first insert Eq. (69) into $`U(𝐯_s)`$ defined in Eqs. (A36), (A37), (A29), (A32), and (A40) to find that

$`U_1^{}(𝐯_s)=8U_2^{}(𝐯_s)=32U_3^{}(𝐯_s)={\displaystyle \frac{2ϵ_F^2\mathrm{sin}4\theta }{m^2r^4}},`$ (72)
$`U_1^{\prime \prime }(𝐯_s)={\displaystyle \frac{1}{96}}U_2^{\prime \prime }(𝐯_s)={\displaystyle \frac{ϵ_F^2\mathrm{cos}4\theta }{4m^2r^4}}.`$ (73)
Equations (60), (72) and (73) together with the asymptotic behaviors of $`L_{xy}`$ in Tabel II imply that in the single vortex case,

$`L_{xy}^{}Y^{}(T)\left({\displaystyle \frac{\xi }{r}}\right)^4\mathrm{sin}4\theta ,L_{xy}^{\prime \prime }Y^{\prime \prime }(T)\left({\displaystyle \frac{\xi }{r}}\right)^4\mathrm{cos}4\theta ,`$ (74)
where $`Y^{}7.4`$ in the GL regime, and $`7.1\mathrm{\Delta }_0(T)/T`$ for $`T\mathrm{\Delta }_0`$, while $`Y^{\prime \prime }9.8`$ in the GL regime, and $`25.9\mathrm{\Delta }_0(T)/T`$ for $`T\mathrm{\Delta }_0`$. Equations (60), (74) and the asymptotic behaviors of $`\eta _{xy}`$ in Tabel II yield

$`\overline{𝒟}_{xy}^{\mathrm{GL}}(𝐫)3.7\left[c_{xy}{\displaystyle \frac{\mathrm{\Delta }_0^2(T)}{8(\pi T)^2}}+{\displaystyle \frac{ϵ_Fmv_s^2(𝐫)}{(\pi T)^2}}\right]^1\left({\displaystyle \frac{\xi }{r}}\right)^4\mathrm{sin}4\theta ,`$ (75)
$`\overline{𝒟}_{xy}^{\prime \prime \mathrm{GL}}(𝐫)4.9\left[c_{xy}{\displaystyle \frac{3}{8}}{\displaystyle \frac{\mathrm{\Delta }_0^2(T)}{(\pi T)^2}}+{\displaystyle \frac{ϵ_Fmv_s^2(𝐫)}{(\pi T)^2}}\right]^1\left({\displaystyle \frac{\xi }{r}}\right)^4\mathrm{cos}4\theta ,`$ (76)
and hence for $`c_{xy}1`$

$`\overline{𝒟}_{xy}^{\mathrm{GL}}(𝐫)0.61ic_{xy}^1\left({\displaystyle \frac{\xi }{r}}\right)^4(e^{i4\theta }+7e^{i4\theta }).`$ (77)
This result coincides with that of Ref. .

At low $`T`$, we obtain

$`\overline{𝒟}_{xy}^{\mathrm{low}T}(𝐫)`$ $``$ $`3.5\left[c_{xy}+{\displaystyle \frac{ϵ_Fmv_s^2(𝐫)}{4\mathrm{\Delta }_0T}}\right]^1{\displaystyle \frac{\mathrm{\Delta }_0}{T}}\left({\displaystyle \frac{\xi }{r}}\right)^4\mathrm{sin}4\theta ,`$ (78)
$`\overline{𝒟}_{xy}^{\prime \prime \mathrm{low}T}(𝐫)`$ $``$ $`12.9\left[c_{xy}{\displaystyle \frac{1}{2}}+2(\mathrm{ln2}){\displaystyle \frac{T}{\mathrm{\Delta }_0}}+{\displaystyle \frac{ϵ_Fmv_s^2(𝐫)}{4\mathrm{\Delta }_0T}}\right]^1{\displaystyle \frac{\mathrm{\Delta }_0}{T}}\left({\displaystyle \frac{\xi }{r}}\right)^4\mathrm{cos}4\theta ,`$ (79)
leading, for $`c_{xy}1`$, to

$`\overline{𝒟}_{xy}^{\mathrm{lowT}}(𝐫)4.69ic_{xy}^1\left({\displaystyle \frac{\mathrm{\Delta }_0}{T}}\right)\left({\displaystyle \frac{\xi }{r}}\right)^4(e^{i4\theta }+1.76e^{i4\theta }).`$ (80)
In Fig. 6, we show the winding of $`d_{xy}`$ component and the magnitude of the normalized $`\overline{𝒟}_{xy}`$ at high $`T`$ (GL) and low $`T`$ respectively. The eight-fold symmetry is more obvious at low $`T`$.

Since many current theories postulate a homogeneous $`d_{x^2y^2}+id_{xy}`$ state without justification, it is interesting to consider the size of the spatially varying part of $`𝒟_{xy}`$ relative to its homogeneous component. The smallest relevant value of the ratio $`𝒟_{xy}^{F_4}/𝒟_{xy}^{OZ}`$, where the superscripts $`F_4`$ and $`OZ`$ indicate the relevant driving terms in Eqs. (60) and (64), respectively, is attained at $`r=R_H`$ in the physically relevant regime $`R_H<r^{}`$. The value is $`(E_F/\mathrm{\Delta }_0)^2(H/H_{c2})`$ in the GL regime, and $`(E_F/\mathrm{\Delta }_0)^2(H/H_{c2})(\mathrm{\Delta }_0/T)`$ at low temperatures; a simple estimate then shows that the spatially fluctuating component is always at least an order of magnitude larger than the homogeneous component at experimentally relevant temperatures and fields.

## V Is there magnetic field-induced pseudo-phase transition?

The results presented in Sections III and IV suggest that there always exist field-induced $`s`$ and $`d_{xy}`$ components in a parent $`d_{x^2y^2}`$-wave superconducting state where the electronic interactions in the $`s`$ and $`d_{xy}`$ channels are nonzero. These subdominant components are spatially inhomogeneous and, in a general case of small $`V_s`$ and $`V_{xy}`$ compared with $`V_d`$, are negligibly small as far as any bulk physical quantity is concerned. However, there is also a special situation in which either $`V_s`$ or $`V_{xy}`$ is nearly degenerate with $`V_d`$, leading to a small $`c_s`$ or $`c_{xy}`$. In this case, the denominator(s) of $`𝒟_i,i=s,d_{xy}`$ may vanish at some critical temperature $`T_{ci}^{}`$, a singularity marking a second pseudo-phase transition into a $`d_{x^2y^2}+𝒟_i`$ state with $`𝒟_i`$ a homogeneous bulk quantity. Investigation of such a possible pseudo-phase transition is particularly interesting in association with the experimentally observed thermal-conductivity plateau as mentioned in the Introduction. Since in high-$`T_c`$ cuprates no such phase transition has been reported to be found in the absence of magnetic field, we focus on the question whether there can be magnetic field-driven phase transition. This is equivalent to searching for a nonzero $`T_{ci}^{}`$ at a finite field which vanishes at zero field.

We first study the $`s`$ component. From Eqs. (51), (52), (54), and (55), we see that the denominator of $`𝒟_s^{}`$ can never be zero, and the singularity may occur in $`𝒟_s^{\prime \prime }`$. In the GL regime, the critical temperatures are

$`T_{cs}^{}(𝐯_s)={\displaystyle \frac{1}{\pi }}\sqrt{c_s^1[1/42ϵ_Fmv_s^2/\mathrm{\Delta }_0^2(T_{cs}^{})]}\mathrm{\Delta }_0(T_{cs}^{}).`$ (81)
It follows that $`T_{cs}^{}(𝐯_s)<T_{cs}^{}(0)`$, implying there is no instability at nonzero field leading to a homogeneous subdominant $`𝒟_i`$. At low $`T`$, $`T_{cs}^{}`$ becomes

$`T_{cs}^{}(𝐯_s)={\displaystyle \frac{1}{4(\mathrm{ln}\mathrm{\hspace{0.17em}2})}}\left[12c_s{\displaystyle \frac{ϵ_Fmv_s^2}{2T_{cs}^{}\mathrm{\Delta }_0}}\right]\mathrm{\Delta }_0.`$ (82)
Again, $`T_{cs}^{}(𝐯_s)<T_{cs}^{}(0)`$ is found, meaning that the magnetic field does not favor such a second-order phase transition from a $`d_{x^2y^2}`$ to $`d_{x^2y^2}+is`$ state. Accounting for nonlinear effects does not affect this conclusion.

A similar analysis can be made in the $`d_{xy}`$ situation. From the asymptotic behaviors of $`\eta _{xy}^{}`$ and $`\eta _{xy}^{\prime \prime }`$ listed in Table II, we find that there is also a singularity in $`𝒟_{xy}^{}`$ in the GL regime, but its corresponding critical temperatures are smaller than that in $`𝒟_{xy}^{\prime \prime }`$. Except for this point, we reach the same conclusion as the $`s`$ case: No magnetic field-induced phase transition is found.

At this stage, we may examine the prospects of a possible low-$`T`$ field-induced first-order phase transition into a time-reversal symmetry breaking state $`d_{x^2y^2}+id_{xy}`$ as suggested by Laughlin . The free energy functional of the $`d_{xy}`$ pairing order parameter assumed by Laughlin takes the form

$`{\displaystyle \frac{F}{L^2}}={\displaystyle \frac{1}{6\pi }}{\displaystyle \frac{(𝒟_{xy}^{\prime \prime })^3}{\nu ^2}}{\displaystyle \frac{1}{\pi }}{\displaystyle \frac{eB}{c}}𝒟_{xy}^{\prime \prime }\mathrm{tanh}^2\left({\displaystyle \frac{𝒟_{xy}^{\prime \prime }}{2T}}\right){\displaystyle \frac{4}{\pi }}{\displaystyle \frac{T^3}{\nu ^2}}\left\{{\displaystyle \frac{(𝒟_{xy}^{\prime \prime })^2}{2T^2}}\mathrm{ln}[1+e^{𝒟_{xy}^{\prime \prime }/T}]+{\displaystyle _{𝒟_{xy}^{\prime \prime }/T}^{\mathrm{}}}\mathrm{ln}[1+e^x]x𝑑x\right\},`$ (83)
where $`\nu =\sqrt{2\mathrm{\Delta }_0/m}`$. The special $`T`$ dependences of the second and third terms lead to a weak first-order phase transition at $`T0.52\nu \sqrt{2eB/c}`$. We may check the justification of this free energy based on the present theory. The first and third terms in the right-hand side of Eq. (83) are understood to come from the local free energy $`F^{(loc)}`$ in Eq. (A5). The second one, as the coupling of the unusual magnetization to the field, is naturally related to the orbital Zeeman term $`F^{OZ}`$ in Eq. (58) with $`M(𝐫)`$ shown in Eq. (B14). But the assumed temperature dependence crucial for the first-order phase transition, is inconsistent with that of $`F^{OZ}`$, which we have shown to be a weak function of $`T`$ for $`T\mathrm{\Delta }_0`$ (see Fig. 1) and a linear function of $`𝒟_{xy}^{\prime \prime }`$ in the limit $`𝒟_{xy}^{\prime \prime }0`$. Therefore, the microscopic calculation in the present work does not confirm Laughlin’s free energy functional which formed the basis of the first-order phase transition found in his work. In the present theory, a small homogenous $`𝒟_{xy}^{\prime \prime }`$ is found for any magnetic field and temperature in the superconducting state. It is swamped by a larger (but still much smaller than $`\mathrm{\Delta }_0`$) $`𝒟_{xy}^{\prime \prime }`$ spatially fluctuating component over almost all of the vortex lattice for physically relevant fields.

## VI Conclusions

In this paper, we have formulated a perturbation theory to investigate the magnetic-field induced subdominant order parameters of a clean $`d`$-wave superconductor in the presence of the gauge-invariant spatially varying supercurrent field in the mixed state. With the assumption of slowly spatially varying supercurrents and their induced $`s`$ and $`d_{xy}`$ components in the bulk sample, we are able to derive the free energy as power series in the Doppler energy shift, the derivatives of the supercurrents, and the subdominant components. The free energy is valid from $`T_c`$ down to very low temperatures, enabling us to compare the resulting $`s`$ and $`d_{xy}`$ components at low $`T`$ with the existing results in the GL regime. To leading order, the real and imaginary parts of the $`s`$ component, driven by local $`m(v_{sx}^2v_{sy}^2)`$ and the derivative of the supercurrent $`_xv_{sx}`$, respectively, were shown to have very different temperature dependences. In the GL regime, both the real and imaginary parts are in the same order. But at low $`T`$, the real part acquires an extra small prefactor $`T/\mathrm{\Delta }_0`$ and the resulting $`s`$ winding happens within a very small region near antinodal directions, leaving a rigid $`d_{x^2y^2}+is`$ state over most of the vortex. It is important to note, however, that this structure does not imply a gap in the quasiparticle spectrum, due to the small size of the $`s`$ component compared with the large Doppler shifts near the core.

To leading order, the $`d_{xy}`$ component is driven by two terms of different symmetries competing over different parts of the vortex lattice. The first is the orbital Zeeman term $`F^{OZ}`$ arising from the coupling of the spontaneous magnetization to the magnetic field. Its significance is limited by its small magnitude due to the particle-hole asymmetric effects. The second driving term $`F_4`$ scales as $`(\xi /r)^4`$ for $`\xi rR_H`$. The crossover scale $`r^{}`$ divides the bulk region into inner and outer regions dominated by the two distinct physical effects. In the inner region, $`F_4`$ determines the vortex structure, and the $`d_{xy}`$ component has eight-fold symmetry and its relative phase winds 4 times that of the $`d_{x^2y^2}`$ component. In the outer region, $`F^{OZ}`$ is more important, leading to a rigid $`d_{x^2y^2}+id_{xy}`$ superconducting state, where the $`d_{xy}`$ component is spatially nearly homogeneous and a weak function of temperature. We have shown that the crossover scale is of order $`r^{}=\xi (R_H/\lambda _F)^{1/2}`$ in the GL regime, and becomes $`(\mathrm{\Delta }_0/T)^{1/4}`$ times larger at low $`T`$. Our best estimate for the high-$`T_c`$ cuprates suggests that $`r^{}`$ is greater than the intervortex separation $`R_H`$ for fields above $`0.01H_{c2}`$, so that it appears that vortex lattice structure in fields of order Tesla is governed by $`F_4`$ and the orbital Zeeman effect is irrelevant. The relative phase between the $`d_{x^2y^2}`$ and $`d_{xy}`$ components, as well as the magnitude of the induced $`d_{xy}`$, are therefore strongly space dependent everywhere in the sample.

Our results have implications for several scenarios which have been proposed to create a state without quasiparticle excitations at low temperatures via the creation of a finite out-of-phase subdominant pair component. No such bulk state is found for generic values of the pair potentials $`V_s`$ and $`V_d`$. It remains possible that induced core $`d_{xy}`$ patches overlap with increasing field, as proposed by Ramakrishnan, leading to a gapped state at high fields beyond the scope of our analysis. However, analyzing the interplay between order parameter suppression $`𝒟_d(𝐫)`$ around the vortex cores and $`𝐯_s(𝐫)`$, we found it to be a very small effect on the induced $`d_{xy}`$ component; it therefore seems unlikely that the field scale where such a transition may occur can be significantly less than $`H_{c2}`$.

Finally, we searched for a possible second order magnetic field-induced phase transition for special values of the coupling constants, and reached a negative conclusion. No such phase transition into a bulk $`d_{x^2y^2}+s`$ or a $`d_{x^2y^2}+d_{xy}`$ state is found unless the transition has already taken place in the absence of the field, which is apparently not the case in the high-$`T_c`$ cuprates. Examining the Laughlin free energy driving a first-order phase transition,  we found that the crucial field-dependent term has an assumed temperature dependence inconsistent with the BCS theory; thus this phase transition picture is not supported in the present work, consistent with numerical results of Yasui and Kita.

The apparent phase transition observed by Krishana et al.  has not been clearly reproduced by other groups, and the possibility exists that the effect is due to inhomogeneously trapped flux. It is still interesting, however, to ask what kinds of intrinsic phase transitions might be possible in a $`d`$-wave superconductor. We have shown that it is unlikely that any bulk phase transition can be induced by a magnetic field, at least in the low-field regime where our approach is valid. Since our model neglects the vortex core regions, it is conceivable that vortex core transitions such as those observed in the $`{}_{}{}^{3}He`$ system might still be relevant. However, since the number of bound states in the core region is small and possibly zero for the present case, it seems unlikely that such a transition would have an important effect on the quasiparticles responsible for heat transport at low $`T`$. A final possibility, currently under investigation within the present framework, is that transitions occur in the vortex lattice structure as a function of field.

## ACKNOWLEDGMENTS

The authors are grateful to N. Andrei, M. Fogelström, M. Franz, C. D. Gong, T. Kita, T. Kopp, S. Sachdev, S. H. Simon, and Y.-J. Wang for helpful communications. Partial support was provided by the A. v. Humboldt Foundation, NSF Grants No. DMR-9974396 and INT-9815833, and “Graduiertenkolleg Anwendungen der Supraleitung” of the Deutsche Forschungsgemeinschaft.

## A Derivation of Free Energy

In this Appendix we derive the free energy $`\stackrel{~}{F}`$ in Eqs. (33-38) which is valid for space region where $`𝐯_s`$ varies slowly. We begin with rewriting Tr$`\widehat{p}^m`$ in Eq. (20) according to Eqs. (30) and (32),

$`\mathrm{Tr}\widehat{p}^m(\mathrm{Tr}\widehat{p}^m)^{(\mathrm{loc})}+(\mathrm{Tr}\widehat{p}^m)^{(\mathrm{der})},`$ (A1)
where $`(\mathrm{Tr}\widehat{p}^m)^{(\mathrm{loc})}`$ corresponds to $`𝐯_s(𝐫_j)`$ and $`𝒟(𝐫_j)`$ for all $`j=1,\mathrm{},m1`$ taking the local values $`𝐯_s(𝐫)`$ and $`𝒟(𝐫)`$, respectively, and becomes

$`(\mathrm{Tr}\widehat{p}^m)^{(\mathrm{loc})}={\displaystyle \frac{d^2𝐤}{(2\pi )^2}d^2𝐫\mathrm{Tr}\left[\widehat{p}_{𝐤,𝐤}(𝐫)\right]^m},`$ (A2)
and $`(\mathrm{Tr}\widehat{p}^m)^{(\mathrm{der})}`$ contains derivatives of $`𝐯_s(𝐫_j)`$ and/or $`𝒟(𝐫_j)`$. Inserting Eq. (A1) into Eq. (20) leads to

$`FF^{(\mathrm{loc})}+F^{(\mathrm{der})},`$ (A3)
where $`F^{(\mathrm{loc})}`$ is just the local free energy obtained in the semiclassical approximation,

$`F^{(\mathrm{loc})}`$ $`=`$ $`F_0T\mathrm{Tr}\mathrm{ln}\widehat{M}_0+T{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{m}}(\mathrm{Tr}\widehat{p}^m)^{(\mathrm{loc})}=F_0{\displaystyle d^2𝐫T\underset{n}{}\frac{d^2𝐤}{(2\pi )^2}\mathrm{ln}[\stackrel{~}{W}_{n𝐤}(𝐫)\eta _𝐤(𝐫)]},`$ (A4)
$`=`$ $`F_0{\displaystyle d^2𝐫\frac{d^2𝐤}{(2\pi )^2}\sqrt{\stackrel{~}{E}_𝐤^2(𝐫)+\eta _𝐤(𝐫)}}{\displaystyle d^2𝐫\frac{d^2𝐤}{(2\pi )^2}\underset{l=\pm 1}{}\mathrm{ln}\left\{1+e^{[l𝐯_s(𝐫)𝐤_F+\sqrt{\stackrel{~}{E}_𝐤^2(𝐫)+\eta _𝐤(𝐫)}]/T}\right\}},`$ (A5)
with

$`\stackrel{~}{W}_{n𝐤}(𝐫)`$ $`=`$ $`[i\omega _n+𝐯_s(𝐫)𝐤_F]^2+\stackrel{~}{E}_𝐤^2(𝐫),`$ (A6)
$`\stackrel{~}{E}_𝐤(𝐫)`$ $`=`$ $`\sqrt{ϵ_𝐤^2+|\overline{\mathrm{\Delta }}(𝐫)|^2\mathrm{\Phi }_{d𝐤}^2}`$ (A7)
$`\eta _𝐤(𝐫)`$ $``$ $`{\displaystyle \underset{i=s,d_{xy}}{}}\{[\overline{\mathrm{\Delta }}(𝐫)𝒟_i^{}(𝐫)+h.c.]\mathrm{\Phi }_{d𝐤}\mathrm{\Phi }_{i𝐤}+|𝒟_i(𝐫)|^2\mathrm{\Phi }_{i𝐤}^2\}.`$ (A8)
Expanding $`F^{(\mathrm{loc})}`$ in power series in $`𝒟_i`$ will give linear-$`\overline{𝒟}_i^{}`$ term as well as quadratic of $`𝒟_i`$ terms. The prefactors are not universal in the whole temperature regimes, as will be discussed and shown below. $`F^{(\mathrm{der})}`$ in Eq. (A3) reads

$`F^{(\mathrm{der})}=T{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{m}}(\mathrm{Tr}\widehat{p}^m)^{(\mathrm{der})}={\displaystyle \underset{j}{}}F^{(j)},`$ (A9)
where we expanded $`(\mathrm{Tr}\widehat{p}^m)^{(\mathrm{der})}`$ as power series in the $`j`$th derivatives of $`𝐯_s`$ or $`𝒟_i`$ with respect to $`𝐫`$: $`(\mathrm{Tr}\widehat{p}^m)^{(der)}=_j(\mathrm{Tr}\widehat{p}^m)^{(j)}`$. $`F^{(\mathrm{der})}`$ reflects nonlocal couplings of the subdominant order parameters to the supercurrent fields. We examine the formal leading term, $`F^{(1)}`$ which includs $`_𝐫v_s(𝐫)`$ or $`_𝐫𝒟_i(𝐫)`$, in Eq. (A9),

$`F^{(1)}={\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{m}}(\mathrm{Tr}\widehat{p}^m)^{(1)}={\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{m}}(\rho _{m1}+\rho _{m2}+\rho _{m3}),`$ (A10)
$`\rho _{m1}=i{\displaystyle \underset{\nu =0}{\overset{m2}{}}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}d^2𝐫\mathrm{Tr}\left\{[_𝐫\widehat{p}_{𝐤,𝐤}(𝐫)]\widehat{p}_{𝐤,𝐤}^{m2\nu }(𝐫)[_𝐪\widehat{p}_{𝐪,𝐤}(𝐫)]_{𝐪=𝐤}\widehat{p}_{𝐤,𝐤}^\nu (𝐫)\right\}},`$ (A11)
$`\rho _{m2}=i{\displaystyle \underset{\nu =0}{\overset{m2}{}}}\nu {\displaystyle \frac{d^2𝐤}{(2\pi )^2}d^2𝐫\mathrm{Tr}\left\{[_𝐫\widehat{p}_{𝐤,𝐤}(𝐫)]\widehat{p}_{𝐤,𝐤}^{m2\nu }(𝐫)[_𝐤\widehat{p}_{𝐪,𝐤}(𝐫)]\widehat{p}_{𝐤,𝐤}^\nu (𝐫)\right\}},`$ (A12)
$`\rho _{m3}=i(m1){\displaystyle \frac{d^2𝐤}{(2\pi )^2}d^2𝐫\mathrm{Tr}\left\{[(_𝐫_𝐪)\widehat{p}_{𝐤,𝐪}(𝐫)]_{𝐪=𝐤}\widehat{p}_{𝐤,𝐤}^{m1}(𝐫)\right\}}.`$ (A13)
$`F^{(1)}`$ can be expanded as power series in $`𝒟_i`$. Note in terms including $`_𝐫𝒟_i(𝐫)`$ the derivative can be transfered to that of $`𝐯_s`$ by partial integral. The resulting linear-$`𝒟_i`$ term in $`F^{(1)}`$ is $`F_𝒟^{(1)}=_{i=d,s,d_{xy}}F_i^{(1)}`$, where

$`F_i^{(1)}`$ $``$ $`i{\displaystyle d^2𝐫T\underset{n}{}\frac{d^2𝐤}{(2\pi )^2}\mathrm{\Phi }_{i𝐤}\{_𝐫[𝐯_s(𝐫)𝐤_F]\}\underset{m=2}{\overset{\mathrm{}}{}}[𝐯_s(𝐫)𝐤_F]^{m2}\underset{\nu =1}{\overset{m1}{}}\nu \mathrm{Tr}\left\{(_𝐤\widehat{g}_𝐤)\widehat{g}_𝐤^\nu \widehat{𝒟}_i(𝐫)\widehat{g}_𝐤^{m\nu 1}\right\}}`$ (A14)
$`=`$ $`{\displaystyle d^2𝐫𝒟_i^{\prime \prime }(𝐫)T\underset{n}{}\frac{d^2𝐤}{(2\pi )^2}\mathrm{\Phi }_{i𝐤}2\mathrm{\Delta }_𝐤[_𝐤ϵ_𝐤_𝐫][𝐯_s(𝐫)𝐤_F]Z_{n𝐤}(𝐯_s)}`$ (A16)
$`{\displaystyle d^2𝐫𝒟_i^{\prime \prime }(𝐫)T\underset{n}{}\frac{d^2𝐤}{(2\pi )^2}\mathrm{\Phi }_{i𝐤}2ϵ_𝐤[_𝐤\mathrm{\Delta }_𝐤_𝐫][𝐯_s(𝐫)𝐤_F]Z_{n𝐤}(𝐯_s)},`$
with $`\widehat{𝒟}_i(𝐫)=\left(\begin{array}{cc}0\hfill & 𝒟_i(𝐫)\hfill \\ 𝒟_i^{}(𝐫)\hfill & 0\hfill \end{array}\right)`$ and

$`Z_{n𝐤}(𝐯_s)={\displaystyle \frac{\left\{[𝐯_s(𝐫)k_F]^2W_{n𝐤}\right\}^24\omega _n^2[𝐯_s(𝐫)k_F]^2}{\left\{[(𝐯_s(𝐫)k_F)^2W_{n𝐤}]^2+4\omega _n^2[𝐯_s(𝐫)k_F]^2\right\}^2}}.`$ (A17)
It is easy to see that $`F_d^{(1)}`$ is negligibly small, so $`F_{\overline{\mathrm{\Delta }}}`$ in Eq. (34) is simply obtained from $`F^{(loc)}`$ in Eq. (A4) by setting $`𝒟_s,𝒟_{xy}=0`$. We argue that $`F^{(loc)}`$ and $`F_s^{(1)}`$ in Eqs. (A4) and (A16) are sufficient for studying the $`s`$ component up to leading orders. Expanding $`F^{(loc)}`$ as power series in $`𝒟_s`$ and $`v_sk_F`$, we will see that the driving term for $`𝒟_s^{}`$ is in leading order of $`m(v_{sx}^2v_{sy}^2)`$, which, in the spatial regime of interest $`\xi <r<R_H`$, scales as $`1/r^2`$. $`F_s^{(1)}`$ gives a driving term for $`𝒟_s^{\prime \prime }`$ scaling as $`_xv_{sx}1/r^2`$. Clearly, terms contained in $`_{l>1}F^{(l)}`$ are of higher orders of $`v_s`$ and derivatives of $`v_s`$. We first show results for $`s`$ component. The $`d_{xy}`$ situation is more complicated and will be discussed later.

### 1 Scaling expressions for $`L_s(T,𝐯_s)`$ and $`\eta _s(T,𝐯_s)`$ at $`T>𝐯_s𝐤_F`$

In this temperature regime, we can expand $`\stackrel{~}{W}_{n𝐤}(𝐫)`$ in Eqs. (A4) and (A16) as power series in $`v_sk_F`$, $`𝒟_i^{}(𝐫)`$, and $`𝒟_i^{\prime \prime }(𝐫)`$ to find explicit expressions for $`L_s(T,𝐯_s)`$ and $`\eta _s(T,𝐯_s)`$ in Eq. (35),

$`L_s^{}(T,𝐯_s)=N_0^1T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\frac{2\mathrm{cos}^22\phi (4\omega _n^2W_{n𝐤})}{W_{n𝐤}^3}ϵ_Fm[v_{sx}^2v_{sy}^2]}=2{\displaystyle \frac{ϵ_Fm[v_{sx}^2v_{sy}^2]}{T^2}}h_{1s}\left(d\right),`$ (A18)
$`L_s^{\prime \prime }(T,𝐯_s)=T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\frac{4\mathrm{cos}^22\phi }{W_{n𝐤}^2}ϵ_F_xv_{sx}}=4{\displaystyle \frac{ϵ_F_xv_{sx}}{T^2}}h_{2s}\left(d\right),`$ (A19)
$`\eta _s^{}(T,𝐯_s)=\eta _s^{\prime \prime }(T,𝐫)+2d^2h_{2s}(d),`$ (A20)
$`\eta _s^{\prime \prime }(T,𝐯_s)=h_{3s}(d)+{\displaystyle \frac{ϵ_Fmv_s^2(𝐫)}{T^2}}h_{4s}(d),`$ (A21)
where Eq. (40) was used, $`d=\mathrm{\Delta }_0(T)/T`$, and

$`h_{1i}(d)={\displaystyle \frac{8}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle _0^{\pi /4}}𝑑\phi \mathrm{cos}^22\phi \mathrm{\Phi }_{i𝐤}^2{\displaystyle \frac{1}{2\eta }}{\displaystyle \frac{d^2f(\eta )}{d\eta ^2}},h_{2i}(d)={\displaystyle \frac{8}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle _0^{\pi /4}}𝑑\phi \mathrm{cos}^22\phi \mathrm{\Phi }_{i𝐤}^2{\displaystyle \frac{t(\eta )}{4\eta ^3}},`$ (A22)
$`h_{3i}(d)={\displaystyle \frac{8}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle _0^{\pi /4}}𝑑\phi (12\mathrm{cos}^22\phi )\mathrm{\Phi }_{i𝐤}^2{\displaystyle \frac{\mathrm{tanh}(\eta /2)}{2\eta }},h_{4i}(d)={\displaystyle \frac{8}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑x{\displaystyle _0^{\pi /4}}𝑑\phi \mathrm{\Phi }_{i𝐤}^2{\displaystyle \frac{1}{2\eta }}{\displaystyle \frac{d^2f(\eta )}{d\eta ^2}},`$ (A23)
with $`\eta =\sqrt{x^2+d^2\mathrm{cos}^22\phi }`$, $`f(\eta )=1/(e^\eta +1)`$, and $`t(\eta )=12f(\eta )+2\eta df(\eta )/d\eta `$.

### 2  $`L_s(T,𝐯_s)`$ and $`\eta _s(T,𝐯_s)`$ at $`T=0`$

At extremely low $`T`$, the perturbation expansions of the free energy as power series in $`v_s(𝐫)k_F`$ and $`𝒟_i^{}(𝐫)`$ by expanding $`\stackrel{~}{W}_{n𝐤}(𝐫)`$ before integrating over $`𝐤`$ breaks down due to the existence of the zero modes in $`E_𝐤`$. The correct approach is to expand the free energy after doing the integrals over $`𝐤`$. For simplicity, we only show the $`T=0`$ results. It is easy to find that the driving term for $`𝒟_s^{}`$ is negligibly small compared with that for $`𝒟_s^{\prime \prime }`$. We may simply set $`𝒟_s^{}=0`$. The local free energy in Eq. (A4) becomes

$`F_{T=0}^{(loc)}=F_0{\displaystyle d^2𝐫\frac{d^2𝐤}{(2\pi )^2}\sqrt{\stackrel{~}{E}_𝐤^2(𝐫)+\eta _𝐤(𝐫)}}+2{\displaystyle d^2𝐫\frac{d^2𝐤}{(2\pi )^2}Y_𝐤(𝐫)\mathrm{\Theta }(Y_𝐤(𝐫))}`$ (A24)
where $`Y_𝐤(𝐫)=\sqrt{\stackrel{~}{E}_𝐤^2(𝐫)+[𝒟_s^{\prime \prime }(𝐫)]^2}𝐯_s(𝐫)𝐤_F`$ and $`\mathrm{\Theta }(x)`$ is the Heaviside step function. We obtain $`F_{T=0}^{(loc)}=\mathrm{\Delta }F+d^2𝐫[c_s+\eta _s^{\prime \prime }(T=0,𝐯_s)][\overline{𝒟}_s^{\prime \prime }(𝐫)]^2`$, where $`\mathrm{\Delta }F`$ is the local free energy at $`𝒟_s=0`$, and

$`\eta _s^{\prime \prime }(T=0,𝐯_s){\displaystyle \frac{1}{2}}+S_\theta {\displaystyle \frac{v_s(𝐫)k_F}{2\mathrm{\Delta }_0}},`$ (A25)
with $`S_\theta =_{l=\pm 1}\left|\mathrm{cos}(\theta +l\frac{\pi }{4})\right|`$.

From $`F_s^{(1)}`$ in Eq. (A16) we get the driving term for $`𝒟_s^{\prime \prime }`$, from which we find

$`L_s^{\prime \prime }(T=0,𝐯_s){\displaystyle \frac{2ϵ_Fm[v_{sx}^2v_{sy}^2]}{\mathrm{\Delta }_0^2}}\left[1S_\theta {\displaystyle \frac{v_s(𝐫)k_F}{\mathrm{\Delta }_0}}\right].`$ (A26)
### 3 Free energy with respect to $`d_{xy}`$ component

We first analyze the driving terms at $`T>𝐯_s𝐤_F`$. The leading-order linear-$`𝒟_{xy}`$ term from $`F^{(loc)}`$ in Eq. (A4) turns out to be

$`F_4^{(loc)}{\displaystyle d^2𝐫\mathrm{\Delta }_0^2(T)N_0Q_3^{}(T)U_3^{}(𝐯_s)\overline{𝒟}_{xy}^{}},`$ (A27)
where

$`Q_3^{}(T)`$ $`=`$ $`N_0^1T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\mathrm{sin}^24\phi \left\{\frac{5}{2W_{n𝐤}^3}+\frac{10E_𝐤^2}{W_{n𝐤}^4}\frac{8E_𝐤^4}{W_{n𝐤}^5}\right\}},`$ (A28)
$`U_3^{}(𝐯_s)`$ $`=`$ $`4ϵ_F^2m^2\left[v_{sx}^2(𝐫)v_{sy}^2(𝐫)\right]v_{sx}(𝐫)v_{sy}(𝐫).`$ (A29)
For $`\xi <r<R_H`$, $`F_4^{(loc)}`$ is of order $`1/r^4`$. As for $`F_{xy}^{(1)}`$ in Eq. (A16), the two terms in the right-hand side are qualitatively different. We leave the discussion on the second term for a bit later. The first term is nonzero in the leading order of $`(mv_s^2)_xv_{sx}`$ which scales as $`1/r^4`$ too. For $`T>v_sk_F`$, this term is

$`F_{xy}^{(11)}{\displaystyle d^2𝐫\mathrm{\Delta }_0^2(T)N_0Q_1^{\prime \prime }(T)U_1^{\prime \prime }(𝐯_s)\overline{𝒟}_{xy}^{\prime \prime }},`$ (A30)
where

$`Q_1^{\prime \prime }(T)`$ $`=`$ $`N_0^1T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\mathrm{sin}^24\phi \left\{\frac{5}{W_{n𝐤}^3}+\frac{6E_𝐤^2}{W_{n𝐤}^4}\right\}},`$ (A31)
$`U_1^{\prime \prime }(𝐯_s)`$ $`=`$ $`ϵ_F^2\left\{m\left[v_{sx}^2(𝐫)v_{sy}^2(𝐫)\right]\left[_xv_{sy}(𝐫)+_yv_{sx}(𝐫)\right]+4mv_{sx}(𝐫)v_{sy}(𝐫)_xv_{sx}(𝐫)\right\}.`$ (A32)
A simple analysis shows that there are terms coming from $`F^{(2)}`$ and $`F^{(3)}`$ in Eq. (A9) which are of the same order. After some algebra, we find that the terms of order $`1/r^4`$ in $`F^{(2)}`$ are

$`F_{xy}^{(2)}`$ $``$ $`{\displaystyle d^2𝐫\mathrm{\Delta }_0^2(T)N_0\left\{Q_1^{}(T)U_1^{}(𝐯_s)+Q_2^{}(T)U_2^{}(𝐯_s)\right\}\overline{𝒟}_{xy}^{}},`$ (A33)
where

$`Q_1^{}(T)`$ $`=`$ $`N_0^1T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\mathrm{sin}^24\phi \left\{\frac{3}{W_{n𝐤}^3}\frac{4(E_𝐤^2+3ϵ_𝐤^2)}{W_{n𝐤}^4}+\frac{16ϵ_𝐤^2E_𝐤^2}{W_{n𝐤}^5}\right\}},`$ (A34)
$`Q_2^{}(T)`$ $`=`$ $`N_0^1T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\mathrm{sin}^24\phi \left\{\frac{2}{W_{n𝐤}^3}+\frac{24ϵ_𝐤^2}{W_{n𝐤}^4}\frac{32ϵ_𝐤^2E_𝐤^2}{W_{n𝐤}^5}\right\}},`$ (A35)
$`U_1^{}(𝐯_s)`$ $`=`$ $`ϵ_F^2\left\{v_{sx}(𝐫)\left[_x^2v_{sy}(𝐫)+3_x_yv_{sx}(𝐫)\right]+v_{sy}(𝐫)\left[_y^2v_{sx}(𝐫)+3_x^2v_{sx}(𝐫)\right]\right\},`$ (A36)
$`U_2^{}(𝐯_s)`$ $`=`$ $`ϵ_F^2[_xv_{sx}(𝐫)][_xv_{sy}(𝐫)+_yv_{sx}(𝐫)],`$ (A37)
and in $`F^{(3)}`$ are

$`F_{xy}^{(3)}`$ $``$ $`{\displaystyle d^2𝐫\mathrm{\Delta }_0^2(T)N_0Q_2^{\prime \prime }(T)U_2^{\prime \prime }(𝐯_s)\overline{𝒟}_{xy}^{\prime \prime }},`$ (A38)
where

$`Q_2^{\prime \prime }(T)`$ $`=`$ $`N_0^1T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\mathrm{sin}^24\phi \left\{\frac{1}{2W_{n𝐤}^3}\frac{2ϵ_𝐤^2}{W_{n𝐤}^4}\right\}},`$ (A39)
$`U_2^{\prime \prime }(𝐯_s)`$ $`=`$ $`ϵ_F^2\left[{\displaystyle \frac{1}{m}}_x^3v_{sy}(𝐫){\displaystyle \frac{1}{m}}_y^3v_{sx}(𝐫)+{\displaystyle \frac{6}{m}}_x^2_yv_{sx}(𝐫)\right].`$ (A40)
The sum of $`F_4^{(loc)},F_{xy}^{(11)},F_{xy}^{(2)}`$, and $`F_{xy}^{(3)}`$ leads to $`F_4(𝐫)`$ term in Eq. (38), with

$`L_{xy}^{}(T,𝐯_s)`$ $`=`$ $`Q_1^{}(T)U_1^{}(𝐯_s)+Q_2^{}(T)U_2^{}(𝐯_s)+Q_3^{}(T)U_3^{}(𝐯_s),`$ (A41)
$`L_{xy}^{\prime \prime }(T,𝐯_s)`$ $`=`$ $`Q_1^{\prime \prime }(T)U_1^{\prime \prime }(𝐯_s)+Q_2^{\prime \prime }(T)U_2^{\prime \prime }(𝐯_s).`$ (A42)
The second term on the right-hand side of Eq. (A16) is nonvanishing by noting that $`[_𝐤\mathrm{\Delta }_𝐤_𝐫][𝐯_s(𝐫)k_F]=\mathrm{\Delta }_0\mathrm{sin}2\phi [_𝐫\times 𝐯_s(𝐫)]_z`$ which picks up the magnetic vector potential $`𝐀(𝐫)`$ in $`𝐯_s(𝐫)`$. This term is nothing but $`F^{OZ}`$ in Eq. (37), with

$`Q^{OZ}(T)=N_0^1T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}2ϵ_𝐤\mathrm{sin}^22\phi Z_{n𝐤}(𝐯_s)}N_0^1T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\frac{2ϵ_𝐤\mathrm{sin}^22\phi }{W_{n𝐤}^2}}.`$ (A43)
$`F^{OZ}`$ is very different from $`F_4`$. Discussions about comparative importance of these two terms can be found in Sec. III.C. The integrand of $`Q^{OZ}`$ includes odd power in $`ϵ_𝐤`$, implying nonzero contribution only in a particle-hole asymmetric system. In order to proceed further, we adopt the following density of states near the Fermi surface in the normal state,

$`N(ϵ)N_0+ϵ{\displaystyle \frac{dN(ϵ)}{dϵ}}|_{ϵ=0}N_0(1+\alpha {\displaystyle \frac{ϵ}{ϵ_F}}),`$ (A44)
to take into account both the particle-hole symmetric and asymmetric contributions. Here $`\alpha `$ is of order unity, and $`ϵ`$ is typically order of $`T`$ or $`\mathrm{\Delta }_0`$, implying that the second term on the right-hand side of Eq. (A44) may be negligible when the contribution from the particle-hole symmetric mode does not vanish.

The quadratic free energy terms in Eq. (36) are extracted from $`F^{(loc)}`$, where

$`\eta _{xy}^{}(T,𝐯_s)=\eta _{xy}^{\prime \prime }(T,𝐯_s)+2d^2h_{2dxy}(d),`$ (A45)
$`\eta _{xy}^{\prime \prime }(T,𝐯_s)=h_{3xy}(d)+{\displaystyle \frac{ϵ_Fmv_{sx}^2(𝐫)}{T^2}}h_{4xy}(d),`$ (A46)
with $`h_{idxy}`$ defined in Eqs. (A22) and (A23).

## B Spontaneous magnetization in the $`d_{x^2y^2}+id_{xy}`$ state

To proceed with a general derivation of the spontaneous magnetization, we rewrite the free energy in Eq. (20) as

$`F=F_0T{\displaystyle \underset{n}{}}\mathrm{Tr}\mathrm{ln}\left(\widehat{G}^1\right)T{\displaystyle \underset{n}{}}\mathrm{Tr}\mathrm{ln}\left[1\widehat{G}\left(\begin{array}{cc}\widehat{V}_1& 0\\ 0& \widehat{V}_2\end{array}\right)\right],`$ (B3)
where $`\widehat{G}=(\widehat{M}_0+\widehat{D})^1`$ with $`\widehat{D}`$ the field-induced off-diagonal pairing order parameter matrix. The orbital Zeeman term coming from the linear-$`\widehat{V}`$ term is

$`F^{OZ}`$ $``$ $`T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤d^2𝐤_1}{(2\pi )^4}d^2𝐫d^2𝐫_1e^{i(𝐤_1𝐤)(𝐫_1𝐫)}\left[G_{11𝐤,𝐤_1}(𝐫)+G_{22𝐤,𝐤_1}(𝐫)\right]\left[(𝐫_1𝐫)_𝐫\right][𝐯_s(𝐫)𝐤_F]}`$ (B4)
$``$ $`T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}d^2𝐫\mathrm{lim}_{𝐤_1𝐤}\left\{i_𝐤\left[G_{11𝐤,𝐤_1}(𝐫)+G_{22𝐤,𝐤_1}(𝐫)\right]_𝐫\right\}[𝐯_s(𝐫)𝐤_F]}`$ (B5)
$``$ $`{\displaystyle d^2𝐫\frac{e}{2mc}T\underset{n}{}\frac{d^2𝐤}{(2\pi )^2}\mathrm{lim}_{𝐤_1𝐤}\left\{(i𝐤\times _𝐤)\left[G_{11𝐤,𝐤_1}(𝐫)+G_{22𝐤,𝐤_1}(𝐫)\right]\right\}𝐁(𝐫)}.`$ (B6)
Comparing Eq. (B6) with Eq. (58) we see that the magnetization is just the expectation value of the magnetic moment operator $`(e/2mc)(i𝐤\times _𝐤)`$,

$`𝐌(𝐫)={\displaystyle \frac{e}{2mc}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\mathrm{Tr}\widehat{d}_𝐤^{}(i𝐤\times _𝐤)\widehat{d}_𝐤}={\displaystyle \frac{e}{2mc}}T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\mathrm{lim}_{𝐤_1𝐤}(i𝐤\times _𝐤)\mathrm{Tr}\widehat{G}_{𝐤,𝐤_1}(𝐫)},`$ (B7)
with $`\widehat{d}_𝐤^{}=(c_𝐤,c_𝐤^{})`$. It is easy to check that $`𝐌=0`$ for a pure $`d_{x^2y^2}`$\- ($`d_{xy}`$-) wave superconductor, and also for a $`d_{x^2y^2}+is`$-wave one in the order of linear $`s`$. For a $`d_{x^2y^2}+id_{xy}`$-wave superconductor with $`d_{xy}`$ component perturbatively small compared to the $`d_{x^2y^2}`$ component, we use the Dyson equation,

$`\widehat{G}_{𝐤,𝐤_1}\widehat{g}_𝐤\delta (𝐤𝐤_1)+\widehat{g}_𝐤\left(\begin{array}{cc}0& 𝒟_{xy}(𝐫)\mathrm{\Phi }_{xy\frac{𝐤+𝐤_1}{2}}\\ 𝒟_{xy}^{}(𝐫)\mathrm{\Phi }_{xy\frac{𝐤+𝐤_1}{2}}& 0\end{array}\right)\widehat{g}_{𝐤_1},`$ (B10)
where $`\widehat{g}_𝐤`$ is defined in Eq. (28). Inserting Eq. (B10) into Eq. (B7) yields

$`𝐌(𝐫)`$ $`=`$ $`{\displaystyle \frac{e}{2mc}}iT{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\mathrm{\Phi }_{xy𝐤}\left[𝒟_{xy}(𝐫)g_{4𝐤}+𝒟_{xy}^{}(𝐫)g_{1𝐤}\right](𝐤\times _𝐤)g_{2𝐤}}`$ (B11)
$`=`$ $`{\displaystyle \frac{e}{2mc}}T{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}𝒟_{xy}^{\prime \prime }(𝐫)\frac{4\mathrm{\Delta }_0ϵ_𝐤\mathrm{sin}^22\phi }{W_{n𝐤}^2}},`$ (B12)
which, up to the leading order of $`D_{xy}`$ and $`v_s`$, is consistent with the second term on the right-hand side of Eq. (A16).

A particle-hole asymmetric system is required to obtain a nonvanishing $`𝐌(𝐫)`$. To understand the physics, we note $`𝐌(𝐫)`$ in Eq. (B11) is in fact

$`𝐌(𝐫)`$ $`=`$ $`{\displaystyle \frac{e}{2mc}}{\displaystyle \underset{n,𝐤}{}}2𝒟_{xy}^{\prime \prime }\left\{(u_𝐤^2v_𝐤^2){\displaystyle \frac{1}{i\omega _n+E_𝐤}}(u_𝐤^2v_𝐤^2){\displaystyle \frac{1}{i\omega _nE_𝐤}}\right\}(𝐤\times _𝐤)g_{2𝐤}.`$ (B13)
with $`u_𝐤^2=(1+ϵ_k/E_𝐤)/2`$ and $`v_𝐤^2=(1ϵ_k/E_𝐤)/2`$ measuring the particle and hole populations, respectively. Eq. (B13) can be interpreted as a magnetic moment contributed from number currents of both particles and holes flowing in the opposite directions, which cancel in a particle-hole symmetric system.

The spontaneous magnetization $`M(𝐫)`$ as a full expression of $`\overline{𝒟}_{xy}^{\prime \prime }`$ can also be obtained from Eq. (B7). After some straightforward algebra we get

$`𝐌(𝐫){\displaystyle \frac{eT}{2mc}}{\displaystyle \underset{n}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^2}\frac{4\mathrm{\Delta }_0(T)𝒟_{xy}^{\prime \prime }(𝐫)ϵ_𝐤\mathrm{\Phi }_{xy𝐤}^2}{W_{n𝐤}\{W_{n𝐤}+[𝒟_{xy}^{\prime \prime }(𝐫)]^2\mathrm{\Phi }_{xy𝐤}^2\}}}.`$ (B14)
$`𝐌(𝐫)`$ in Eq. (B14) is important for the direct comparison with Laughlin’s free energy as discussed in Sec. V.