Anomalous Scaling of the Penetration Depth in Nodal Superconductors

Anomalous Scaling of the Penetration Depth in Nodal Superconductors

Abstract

Recent findings of anomalous super-linear scaling of low temperature () penetration depth (PD) in several nodal superconductors near putative quantum critical points suggest that the low temperature PD can be a useful probe of quantum critical fluctuations in a superconductor. On the other hand, cuprates which are poster child nodal superconductors have not shown any such anomalous scaling of PD, despite growing evidence of quantum critical points. Then it is natural to ask when and how can quantum critical fluctuations cause anomalous scaling of PD? Carrying out the renormalization group calculation for the problem of two dimensional superconductors with point nodes, we show that quantum critical fluctuations associated with point group symmetry reduction result in non-universal logarithmic corrections to the -dependence of the PD. The resulting apparent power law depends on the bare velocity anisotropy ratio. We then compare our results to data sets from two distinct nodal superconductors: YBaCuO and CeCoIn. Considering all symmetry-lowering possibilities of the point group of interest, , we find our results to be remarkably consistent with YBaCuO being near vertical nematic QCP, and CeCoIn being near diagonal nematic QCP. Our results motivate search for diagonal nematic fluctuations in CeCoIn.

pacs:

I Introduction

The experimental study of quantum criticality associated with a quantum critical point (QCP) located inside a superconducting phase is challenging due to the dominance of superconducting response. Traditional probes for detecting effects of quantum critical fluctuations, such as transport and specific heat, are frequently overwhelmed by superconductivity. For this reason, recent experiments observing anomalous behavior in the temperature dependence of the penetration depth (PD), which is a direct measure of the superfluid density, near putative QCP’s have raised hope in using this fundamental observable for a superconductor (SC) to probe the effects of quantum critical fluctuations in a superconducting system. In particular several different nodal superconductors including heavy fermions [Ormeno et al., 2002; Chia et al., 2003; Özcan et al., 2003; Hashimoto et al., 2013; Truncik et al., 2013], iron-pnictides [Hashimoto et al., 2013] and organic superconductors [Carrington et al., 1999] have shown super-linear -dependence of low temperature PD (as opposed to -linear dependence expected of nodal superconductors [Hardy et al., 1993]). This unusual temperature dependence in a narrow doping range invites one to invoke quantum criticality.

Proximity to a putative antiferromagnetic QCP in the systems investigated in Ref.[Hashimoto et al., 2013] led Hashimoto et al. to conjecture that the antiferromagnetic quantum critical fluctuations cause the observed superlinear scaling. However, since the antiferromagnetic order parameter field carries a finite momentum that does not nest the nodes, it cannot couple to nodal quasiparticles linearly while preserving crystal momentum. Hence the coupling between nodal quasiparticles and antiferromagnetic quantum critical fluctuations are irrelevant in the renormalization group sense and unlikely to alter the temperature dependence of PD in qualitative manner. Therefore to test the possibility of the quantum critical fluctuation driven anomalous scaling scenario, the search for the candidate quantum critical fluctuations should be broadened.

Figure 1: (Color online) Schematic phase diagram of a quantum phase transition inside the superconducting dome associated with developing a new order represented by an order parameter . Here is the superconducting transition temperature, and is a tuning parameter that drives the system through the quantum critical point at . Near there is a quantum critical region where the penetration depth displays anomalous scaling.

In this paper, we explore other avenues for a QCP inside the superconducting dome (see Fig. 1) altering the temperature dependence of the PD of nodal superconductors. In particular, we will be interested in QPT’s that preserve the nodes [Sato, 2006; Berg et al., 2008]. This narrows the possibilities to time-reversal symmetric point group lowering transitions including different types of nematic order [Fradkin et al., 2010]. These possibilities which consist of QPT’s associated with one-dimensional (Fig. 2(a-c)) and two-dimensional representations (see Fig. 2(d)) of the point group and they exhaust all cases of quantum critical fluctuations that can couple linearly with the nodal quasi-particles through non-derivative coupling.

For the QCP’s associated with the one-dimensional representations(Fig. 2(a-c)), we compute leading temperature dependence of the PD and find logarithmic corrections to the -linear PD expected of nodal superconductors away from QCP’s. For the QCP associated with the two-dimensional representation (see Fig. 2(d)) we only state the form of the action for completeness, without carrying out the full calculation of PD. This decision is based on two reasons: 1) the calculation of temperature dependent PD for such QCP requires an involved calculation due to having two coupling terms that compete with each other, 2) such doubling of nodes, should it happen, will be evidenced more directly via other observables. We then compare our results to the PD measurements of YBaCuO and CeCoIn.

Figure 2: (Color online) Nodal structure of -SC with coexisting orders for a single band system. For multi-band systems like CeCoIn, the nodes on different surfaces will move simultaneously. The dashed (blue) lines at the corner of the Brillouin zone represent the Fermi surface, and the diagonal dashed (black) lines represent where the gap vanishes. The open circles denote the original nodes of -SC, and the filled circles denote the reconstructed nodes with coexisting orders. (a) B axial nematic case; (b) B diagonal nematic case; (c) A case; (d) E case. Note for (a-c), each nodal point is spin degenerate. (a, c) correspond to deformation of the gap (particle-particle channel), and hence are associated with in Nambu space. (b) corresponds to deformation of the Fermi surface (particle-hole channel), and hence . The spin degeneracy is lifted for (d) and the blue and red filled circles have different spin-momentum helicity. The arrows in (d) indicate in-plane spin orientation at the nodal points, and here the ordering is predominantly in the particle-hole channel.

Ii The Model

Although the nodal superconductors studied in Ref. [Ormeno et al., 2002; Chia et al., 2003; Özcan et al., 2003; Hashimoto et al., 2013; Truncik et al., 2013; Carrington et al., 1999] all differ in terms of details of their nodal structure and the origin of nodes, we will focus on models of point nodes in two-dimensional system as a representative example. Furthermore, for concreteness, we consider the nodal superconductor to be a -SC on a two-dimensional square lattice with appropriate number of nodal Dirac fermion species: two per spin for cuprates, six per spin for CeCoIn. Below we first present the single band model with two species of Dirac fermions, and then extend it to the three band model relevant for CeCoIn.

ii.1 Single band model: the case of cuprates

Targeting at cuprates, the QCP associated with the quantum phase transition inside the dome of -SC has been studied extensively using a single band model emplying field theoretical approaches [Vojta et al., 2000; Kim et al., 2008; Huh and Sachdev, 2008; Xu et al., 2008; Fritz and Sachdev, 2009; Wang et al., 2011; Liu et al., 2012; Wang, 2013]. In particular an “infinite anisotropic” fixed point was discovered in Ref. [Kim et al., 2008] associated with nodal nematic QCP, and a full renormalization group (RG) theory for the fixed point was developed in Ref. [Huh and Sachdev, 2008]. Following these references, the effective action for the nodal quasi-particles of this -wave superconductor is obtained by linearizing the Bogoliubov-de Gennes mean field Hamiltonian near the two nodes , and their time-reversal partners. We define the momentum deviation from the nodal points near each node , and the nodal Dirac fermions in terms of the Nambu spinors , where is the antisymmetric Levi-Civita symbol, and are spin indices. The effective action is

(1)

where is the Matsubara frequency, and are the Fermi velocity and gap velocity respectively, and are Pauli matrices in the Nambu space that represent kinetic energy and pairing terms respectively. The momentum has been rotated by to simplify the notation. Appendix I demonstrates how this effective action gives rise to a linear-in- temperature dependence of the PD.

The complete catalogue of time-reversal symmetric point group lowering possibilities must correspond to the non-trivial representations of the point group of interest. For C, these representations are: A, B, B and E, where the first three are one-dimensional and the E representation is two dimensional. Associated with each of these representations, we construct the corresponding order parameter fields whose number of components equals the dimensionality of the representation. For A, B and B, the order parameter actions are all of the form

(2)

though they each couple differently to the quasi-particles. For the E representation case, the order parameter action is similar but with for .

To obtain an intuitive picture of the four types of phase transitions, we plot their effects on the nodes in Fig. 2. For simplicity we are showing the case of single Fermi surface, and it can be easily extended to multiple Fermi surfaces. As can be seen in Fig. 2(a) and (b), the B and B cases are axial and diagonal nematic or orientational orderings respectively. The A case rotates the nodes clockwise or counter-clockwise. Finally, for the E case which breaks inversion symmetry, time reversal symmetry requires a Rashba type spin-orbit coupling that splits the nodes. Mean-field order parameters for shifting nodes along the tangent or normal of the Fermi surface involve bilinears of Nambu spinors combined with or respectively. From this perspective, B(axial nematic) and A channels involve , B (diagonal nematic) channel involves , and the E channel involves both.

Now we turn to the leading couplings between the order parameter fields and the nodal quasiparticles. This coupling is generally of the form

(3)

where we rescaled the bosonic field by the coupling strength, , and used the four-component Nambu spinor . For the scalar field cases, the couplings are diagonal in spin space with for the A case, for the B case and for the B case. We will see that the important distinction between the three cases here is not the signs of the couplings but whether they involve -type coupling that belong to the pairing channel and -type coupling that belong in the particle-hole channel. For the E case, where is a linear combination of , we have , where are the spin Pauli matrices, and denotes the ratio between the couplings in particle-particle and particle-hole channels. With two competing couplings, the calculation is more involved for the E case and we will leave it for a future study.

ii.2 Multiband model: the case of CeCoIn

Now we generalize the above single-band low energy effective model to multi-band systems in order to describe the temperature dependence of penetration depth for CeCoIn at low temperatures. CeCoIn is electronically quasi two-dimensional [Petrovic et al., 2001; Kohori et al., 2001; Sidorov et al., 2002]. Its superconducting gap has point-like nodes, and most likely has symmetry, with point nodes in two dimensions along directions, as shown by the magnetic field-angle dependence of thermal conductivity [Izawa et al., 2001] and specific heat [An et al., 2010], and more recently by STM quasiparticle interference spectroscopy [Allan et al., 2013; Zhou et al., 2013]. Hence the low energy degrees of freedom in the superconducting state of CeCoIn can be modeled by two dimensional Dirac fermions around the nodal points. The electronic structure has symmetry, and hence the corresponding new phases can be classified according to the irreducible representations of as shown above for cuprates.

However, there is one extra complication for CeCoIn as compared to cuprates: there are three sets of nodal points associated with the three different Fermi surfaces 1. Hence the low energy degrees of freedom in the superconducting state involve three sets of nodal Dirac fermions. The fermion action becomes

(4)

with the three sets of velocities , where , labeling different Fermi surfaces. The Dirac fermions couple to the order parameter fluctuations with different coupling strength , with the action

(5)

The bare action of the order parameter field remains the same as in the single band model (Eq. (2)).

order parameter coupling
B, A
B
Table 1: Summary of main results: the representation of the new order parameter, the involved Nambu space Pauli matrix for the coupling, fixed point value of the velocity ratio under RG, and the temperature dependence of penetration depth. Here is the temperature corresponding to the RG scale , which is related to the initial value of the velocity ratio by for -coupling, and for -coupling.

Iii -dependence of the Penetration Depth

The -dependence of PD near the putative QCPs can be determined by the RG equation for the paramagnetic part of the electromagnetic response kernel , with . In order to derive the necessary RG equations, we use the standard expansion which is known to be controlled by the smallness of the velocity ratio (e.g., ) near the infinitely anisotropic fixed points even for realistic value of [Huh and Sachdev, 2008].

iii.1 Single band model

For the single band model, the response kernel satisfies the following simple homogeneity relation (see Appendix II):

(6)

with the energy scale , where is the UV energy scale, are the velocities at scale , and are the corresponding velocities at scale . Eq. (3) states that the scaling of is determined by scaling of the two velocities and . This is a result of the cancellation of the renormalization of the fermion field and fermion-electromagnetic field coupling which prevents renormalization of the coupling between the fermion current and the electromagnetic gauge field. Hence the anomalous dimension for the current-current correlation function vanishes. But by dimensional analysis, , with a scaling function . Then we obtain the RG equation for using the scaling relation Eq.(6) to be:

(7)

This shows that the RG flow of the kernel is entirely determined by the flow of the velocity anisotropy ratio .

The RG flow of the velocity anisotropy ratio was calculated for -coupling (e.g., B and A channel ordering, see Fig2(a), (c)) in [Huh and Sachdev, 2008], and the corresponding RG flow for -coupling (e.g., B channel ordering, see Fig2(b)) can be obtained by a “duality” transformation (exchanging with , and exchanging the functional dependence of their RG coefficients, see supplementary material S2). The fixed points associated with each of these flows are for -coupling, and for -coupling. Using these RG flows we now derive the asymptotic temperature dependence of near the fixed points. Consider first -coupling, where the velocity anisotropy ratio is of the asymptotic form [Huh and Sachdev, 2008], and hence . Carrying out a similar process also for -coupling one obtains

(8)

where the exponential decay comes from the engineering dimension of , and the power law prefactor arising from the coupling of nodal fermions to critical modes. Using the RG equation for temperature, , one then obtains the asymptotic temperature dependence of the kernel near fixed points as shown in Table I. The only free parameter here is the initial value of the velocity ratio.

Figure 3: (Color online) Temperature dependence of the penetration depth (or equivalently the superfluid density with ), fitting experimental results for YBaCuO (a) and CeCoIn (b). The diamonds and squares are the experimental results for YBaCuO from [Hardy et al., 1993] and for CeCoIn from [Truncik et al., 2013] respectively. The experimental values of have been normalized by an overall factor so that . The solid lines are fit to the analytic expressions shown in Table I. The fitting parameters are: (a) for -coupling, as deduced from the experimentally measured velocity anisotropy [Chiao et al., 2000; Norman et al., 1995] with ; (b) for -coupling with , as best fit, and correspondingly for all three bands .

iii.2 Multiband model

The above calculational procedure can be straightforwardly generalized to the multiband models. In the three band model for CeCoIn, there are three coupling constants , , that flow. However, it turns out they flow to a stable fixed point where exhibiting emergent enlarged symmetry (see Appendix II). Therefore the behavior of the system near the fixed point is well captured by a single band model with in the RG equation (7) replaced by a summation of the velocity ratios from different bands:

(9)

where .

The asymptotic form of the velocity ratios remains of the same logarithmic form as in the single band model, since the logarithm comes from integration over the fermion propagator, and the asymptotic form of the fermion propagator does not change with the presence of multiple bands. As a result, the RG equation of the electromagnetic kernel, and consequently the temperature dependence of the penetration depth, are of the same asymptotic form as in the single band model.

iii.3 Comparison with experiments

Our result shows that coupling between nodal quasiparticles and quantum critical fluctuations cause logarithmic corrections to the -linear temperature dependence of a nodal superconductor’s inverse penetration depth rather than a universal power law. Such logarithmic corrections would yield “apparent” power law with sublinear (superlinear) -dependence of the PD depending on coupling in the particle-particle channel (particle-hole channel). The “apparent” exponent will be non-universal and depend on the bare velocity anisotropy ratio.

Now we use the predicted -dependence of the PD (Table I) to fit experimental results for YBaCuO (optimal doping) and for CeCoIn. The temperature dependence of for YBaCuO has been favorably compared with -linear behavior expected for -wave superconductor with nodal quasi-particles, despite growing evidence of quantum critical fluctuations [Daou et al., 2009; Ramshaw et al., 2015]. However fitting the data from Ref. [Hardy et al., 1993] with -coupling (e.g. axial nematic order, B1 case in Fig. 2a), we find the extreme value of velocity anisotropy [Chiao et al., 2000; Norman et al., 1995] implies that the existence of nematic quantum critical fluctuations in YBaCuO cannot be ruled out.

The experimental data in [Truncik et al., 2013] for CeCoIn are best fit by -coupling (e.g., diagonal nematic order, B2 case in Fig. 2b) with moderate velocity anisotropy ratio of for all three bands (see Fig. 3b)2. Thus the observation of superlinear -dependence of PD in CeCoIn is consistent with logarithmic correction to penetration depth due to quantum critical diagonal nematic fluctuations.

Iv Conclusion

In conclusion, we have shown in this paper that inside a two dimensional nodal superconductor with point nodes, quantum critical fluctuations associated with TRS-preserving orders that have nonvanishing amplitude at the nodal points give rise to logarithmic corrections to the temperature dependence of PD (see Table I). In particular, we find the direction of fluctuations in the nodal positions play the key role: fluctuations along the direction normal (tangent) to the underlying Fermi surface cause apparent super- (sub-) linear scaling of the PD as a function of temperature. Our results are qualitatively different from predictions of -linear to -square crossover behavior in PD of different origins, e.g. impurity scattering [Hirschfeld and Goldenfeld, 1993], nonlocal effects [Kosztin and Leggett, 1997], antiferromagnetic fluctuations [Nomoto and Ikeda, 2013].

With its normal state displaying non-Fermi liquid behavior, CeCoIn has long been suggested to be close to a QCP [Petrovic et al., 2001; Kohori et al., 2001; Sidorov et al., 2002; Bianchi et al., 2003]. The QCP is associated with a magnetic phase transition based on empirical evidence: for example, replacing 15-30 In by Cd produces an antiferromagnetic ground state [Pham et al., 2006]. Our analysis of the temperature dependence of PD, and also the studies on zero-temperature PD in [Chowdhury et al., 2013; Levchenko et al., 2013], reveal that the magnetic QCP cannot be the whole story. In fact, our result suggests that CeCoIn may be close to a QCP of diagonal nematic order.

Interestingly, previous theoretical studies on quantum critical diagonal nematic fluctuations [Doh et al., 2006; Kao and Kee, 2007] have shown that continuous quantum phase transition associated with Pomeranchuk instability [Oganesyan et al., 2001] is possible with the diagonal nematic even in the presence of lattice unlike the axial nematic case. Further they also showed that such fluctuations cause non-Fermi liquid behavior.

A direct confirmation of strong nematic fluctuations in CeCoIn would enlarge the class of materials exhibiting point group symmetry breaking. This has so far been established in the cuprates [Lawler et al., 2010; Mesaros et al., 2011], iron-pnictides [Chu et al., 2010; Chuang et al., 2010], and URuSi [Okazaki et al., 2011; Ikeda et al., 2012]. Furthermore, establishment of quantum criticality inside superconducting phase would contribute to the emerging picture of intertwined orders in strongly correlated electron systems [Fradkin et al., 2015; Davis and Lee, 2013]. Hence detection of nematic fluctuations in CeCoIn can offer valuable experimental and theoretical insights. This can be achieved following recent developments in using response to strain to measure nematic fluctuations [Chu et al., 2012; Böhmer et al., 2014; Riggs et al., 2015].

We thank Steve Kivelson for helpful discussions. J.-H.S and E.-A.K were supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Science and Engineering under Award de-sc0010313.

Appendix I: Penetration depth for -wave superconductor

In this section, we derive the -dependence of penetration depth for a -wave superconductor from the action of free Dirac fermions (Eq.(1) in main text). The electromagnetic kernel at zero external frequency can be written as

(10)

with the fermion Green’s functions near the two pairs of nodes and . Here to compare with the result of Kosztin and Leggett [Kosztin and Leggett, 1997], we consider also at finite momenta. The Fermi velocity can be approximated by its value at the corresponding nodal point. The momentum summation is only non-zero when , and . After taking the trace, integrating over momentum, replacing the frequency summation by an integral, and defining , one obtains

(11)

We note that the above result is of the scaling form

(12)

with exponent and dynamical exponent .

For a specular boundary, the temperature dependent part of penetration depth can be determined from

(13)

with [Kosztin and Leggett, 1997]. One thus obtains for a -wave superconductor

(14)

with the crossover temperature . The penetration depth crosses over from -linear at higher temperatures to at lower temperatures, as was shown by Kosztin and Leggett in [Kosztin and Leggett, 1997].

Appendix II: Renormalization of the electromagnetic kernel

We present here the renormalization group calculation of the electromagnetic kernel for the system of nodal Dirac fermions coupled with critical fluctuations. A RG theory of such systems has been developed by Huh and Sachdev in [Huh and Sachdev, 2008] for single band models. Here we extend their approach to multiband systems, and to the calculation of correlation functions. The RG coefficients for the renormalization of the fermion Green’s function have been obtained for -coupling in [Huh and Sachdev, 2008]. The results for -coupling can be obtained in a similar way, which turn out to be related to the case by a “duality” transformation as shown below. The RG coefficient of the current vertex is shown to be the same as the anomalous dimension of the fermion field, which ensures that the coupling between the fermion current and the electromagnetic gauge field is not renormalized.

The electromagnetic kernel is defined as

(15)

where the current is

(16)

To simplify the notation, we keep the band index implicit (e.g. ). Adding electromagnetic field to the system introduces a coupling term to the action and allows us to write as

(17)

with the partiton function in the presence of the electromagnetic field

(18)

Electromagnetic kernel of free fermions from RG

To set the stage, we consider first free Dirac fermions coupled with the external gauge field. The partition function reads

(19)

The RG procedure we follow begins by separating slow and fast modes. For the current this implies

(20)

For the partition function this implies

(21)

Integrating out the fast modes generates a new term in the action that is a functional of ,

(22)

where the leading order term represents the superfluid stiffness

(23)

with

(24)

The partition function then obeys

(25)

Then one obtains from Eq.(17),

(26)

With , where is infinitesimally small, one can expand the above equation to yield

(27)

where we have used . This then reproduces the linear -dependence of the kernel.

Electromagnetic kernel for fermions coupled to critical modes from RG

Now consider coupling the nodal fermions to the critical fluctuations. The partition function can be written as

(28)

Note that since the critical fluctuations have , they do not mix fermions in different bands. The fermion-critical mode coupling now renormalizes the action of the slow modes , and the coupling to gauge field . It also contributes to the superfluid stiffness term via coupling to . We consider these effects one by one.

Consider first the renormalization of the fermionic part of the action. Using the bare action of slow fermions

(29)

with , integrating out fast modes leads to a correction to the effective action

(30)

with the “shell” self energy

(31)

where is the self energy obtained by integrating modes with , and its cutoff dependence can be written as

(32)

The action then reads

(33)

with the renormalized velocities for infinitesimal

(34)
(35)

Then one performs the rescaling,

(36)
(37)
(38)
(39)

To recover the original form of the action, it is required that . The action after rescaling is then

(40)

Consider next the renormalization of the coupling between fermions and the electromagnetic field. The bare coupling is

(41)

which also receives corrections after integrating out fast modes,

(42)

with “shell” vertex correction

(43)

The vertex correction is obtained by integrating modes with , and it can be parameterized as , with the dimensionless part obeying

(44)

The renormalized fermion-electromagnetic field coupling then becomes

(45)

After rescaling, using the relation , which we prove in the following subsection, one obtains

(46)

which is of the same form as the unrenormalized coupling term, in accordance with Ward identity.

Consider finally the superfluid stiffness term generated from integrating out fast modes. We use the standard expansion. Each fermion loop contributes a factor , and each vertex contributes . The leading order diagrams are the ring diagrams in the RPA approximation, which scale as . The next order diagrams are the fermion self-energy and vertex corrections, which scale as . However since the fermion-electromagnetic field coupling is proportional to the identity matrix in Nambu space, and the fermion-critical mode coupling is proportional to or , the two fermion loops at the two ends of the ring diagram, which involve both couplings, vanish upon tracing over the Pauli matrices. The leading order contributions from the critical modes are thus the fermion self-energy and vertex corrections. A resummation of the corresponding Feynman diagrams gives

(47)

with the renormalized Green’s function and vertex

(48)
(49)

Using the relation , one obtains

(50)

The partition function then obeys

(51)

and consequently the kernel

(52)

Expanding the R.H.S. of the above equation for infinitesimal yields the Callan-Symanzik equation

(53)

where

(54)
(55)

The solution for finite is of the form

(56)

But by dimensional analysis,

(57)

Restoring the band index, we have

(58)

where . This equation, together with the RG equations for the velocity ratios and the coupling strengths

(59)
(60)

fully determine the RG flow of the electromagnetic kernel, and hence the temperature dependence of the penetration depth at low temperatures.

RG coefficients

The RG equations for the multiband models obviously involve more variables than that of the single band model and hence are more complicated. However we note that near the fixed point, the behavior of the system simplifies significantly due to the enlarged symmetry at the fixed point. Let us consider for example -coupling with bands, which is relevant for CeCoIn. By solving the RG equations numerically, one can show that the fixed point is at , and . The fixed point has an enlarged symmetry. The asymptotic behavior of such a multiband system near the fixed point is thus the same as that of the single band system, with the number of fermion flavors changed from to . Hence in real calculations, it is sufficient to retain only the RG coefficients in the single band model, which are included below.

For -coupling, the RG coefficients have been calculated explicitly in [Huh and Sachdev, 2008],

(61)
(62)
(63)

with

(64)

The results for -coupling can be obtained by interchanging the two velocities, and with in the above equations, i.e.

(65)
(66)