text:
# Correlations of triggering noise in driven magnetic clusters

## I Introduction

Field driven disordered ferromagnets at low temperatures exhibit Barkhausen (BH) noise, a very important physical phenomenon, which is used for noninvasive characterization technique in commercial alloys. It has been recognized that measured BH noise exhibits scale free behavior when only the external field is varied for all values of disorder and driving fields in real experiments  (a short summary of the experimental data can be found in Ref. ). This is in sharp contrast to some theoretical conclusions  emphasizing fine tuning of the strength of disorder to a single critical value. It has been also attempted to understand the dynamics of domain walls which results in BH avalanches in terms of the dynamics of a sandpile model  (for a recent review of sandpile models see ) and of the models of interface depinning . It should be stressed that in BH noise the presence of disorder plays an important role via pinning of domain walls. It remains to be understood how the motion of domain walls is affected by the pinning.

Different scaling behaviors of BH noise can be attributed to the domain structure, which is related to annealing and type of impurities, and nucleation and coalescence of domains, as well sa varying driving conditions. Annealing the samples in the applied anisotropic stress or in the magnetic field leads to a characteristic structure of extended system-size domains with a $`180^{}`$ domain walls parallel to the anisotropy axis . Numerical simulations with extended domain wall in two dimensions  led to the conclusion that different scaling behavior can be expected in two limiting cases: (a) When disorder is weak BH response is dominated by motion of a single (extended) domain wall; (b) For strong disorder a multidomain structure occurs with many competing domain walls. This picture is in a qualitative agreement with experimental results in stress annealed Fe-B-Si  and Fe-Co-B  alloys. Therefore the two different universality classes can be related to surface and bulk criticality, respectively. Numerical simulations starting from an uniform ground state , on the other hand, do not take into account extended domain walls, and thus correspond to the behavior at high-disorder. In both cases, however, the origin of scaling has not been fully understood.

Recently the exact results of the random-field Ising model on the Bethe lattice  show that the avalanche distributions at a fixed driving field have finite cutoffs. However, an infinite cutoff appears for a range of disorder strengths $`\mathrm{\Delta }<\mathrm{\Delta }_c`$ if a distribution is integrated over the hysteresis loop. The integration thus involves an infinite jump in magnetization (when system size $`L\mathrm{}`$) at a critical value of driving field. On the other hand, for strong disorder $`\mathrm{\Delta }`$ above $`\mathrm{\Delta }_c`$ the avalanche distributions remain exponentially damped for all fields . These results encourage further study of the field-integrated distributions to elucidate the role of driving field in the appearance of the scale-free BH noise.

In this work we study properties of the threshold driving field which triggers Barkhausen avalanches in the multidomain structure, which is generated in the linear part of the ascending branch of hysteresis. When the system is slowly driven by increasing the external magnetic field $`H`$ with time, even by continuous field changes, there exists a threshold field that equals weakest pinning force of a domain wall in the system and thus starts an avalanche. In the numerical experiment we can adjust the driving field updates to the weakest local field in a system (so-called infinitely slow driving), and thus we can examine precisely the properties of the triggering noise. Surprisingly, we find that time series of such field updates are long-range correlated. Both Fourier spectrum of the threshold field fluctuations $`\mathrm{\Delta }H(t)H(t+1)H(t)`$ and distribution of distances of the successive avalanches which that field triggers appear to decay with a power law. The exponents are (weakly) universal in a range of values of disorder where BH avalanches have a power-law distributions. In addition, we find that the distribution of trapping times of domain wall at a given point in space exhibits scaling behavior. We show how the long-range noise correlations and fractal properties of the trapping time distribution are related to the observed scaling behavior of BH avalanches.

The present study is motivated by recently observed pattern formation  and activity correlations along a driven interface  in $`1+1`$ dimensions stacked by random defects. In contrast to these models, here we have a 2-dimensional system with many interacting interfaces, and two separate time scales (compared to extremal driving): slow time scale of field updates, and avalanche propagation time scale between two field updates. Despite of these differences, in both cases the intermittent avalanche-like dynamics occurs, which we believe is essential for the observed scaling behavior.

## II Model and simulations

We consider a simple model with disorder represented by local random fields $`h_x`$, which appears from an original disorder via coarse-graining

$$=\underset{<x,x^{}>}{}J_{x,x^{}}S_xS_x^{}\underset{x}{}(h_x+H)S_x,$$
(1)
where $`x(x_{},x_{})`$ and $`J_{x,x^{}}=1`$ is a constant interaction between nearest-neighbor spins $`S_x=\pm 1`$. A Gaussian distribution of $`h_x`$ is assumed with zero mean and width $`f`$. (Other types of disorder have been also considered ). A domain wall of reversed spins is created along $`<11>`$ direction  on the square lattice rotated by $`\pi /4`$. This model is motivated by the extended domain walls in stress-annealed samples, as discussed in Ref. . Periodic boundaries are applied in the direction of wall and an open boundary at the opposite side of the wall. The system is driven globally—updated value of the external field is applied to all spins in the system. The dynamics consists of spin flips when the local field exceeds zero (see below).

As discussed in Ref. , the fact that motion of the $`11`$ domain wall has no energy threshold at vanishing disorder has several advantages. In particular, the wall depinning occurs along a line $`H_c(f)\kappa f`$ in the $`(f,H)`$ plane, where $`\kappa `$ is a constant (see Fig. 2 in Ref. ). Therefore for the domain wall along $`11`$ direction an infinitesimally small field $`H`$ is sufficient to move and depin the domain wall for small disorder, $`H_c(f)0`$ for $`f0`$. This important property of the model allows us to use much smaller lattice sizes compared to earlier studies , where a large system has to be used in order to find a random field large enough to surmount the energy barrier $`2J`$ in the case of wall along $`10`$ direction, or the nucleation energy $`4J`$ in the case of uniform ground states. On the other hand, by increasing disorder the distance between pinning centers decreases, and thus a critical disorder $`f^{}`$ exists at which depinning is no longer possible. Instead, it becomes energetically easier to nucleate new domains in the bulk. The exact value of $`f^{}`$ depends on the type of distribution , and is still not known exactly. In Ref.  it was found using a finite-size scaling analysis of the averaged interface velocity and avalanche distributions that for the present model $`f^{}=0.62`$ within numerical error bars.

The scaling properties of BH avalanches in low disorder regime were discussed in Ref. . Here we are interested in the high disorder region, where the $`11`$ interface remains pinned at all fields. Nevertheless, the initial state with the $`11`$ interface ensures that the ratio of the characteristic lengths $`\xi _1/L1`$ holds for the applied disorder values (thus the avalanche size cut-off $`s_0L^2`$). (The characteristic length $`\xi _1`$ represents distance between strong pinning centers  for given strength of disorder.) On the contrary, when $`\xi _1/L1`$, the $`11`$ interface moves, meaning that the selected system size $`L`$ is small for given strength of pinning, and thus system is not in the high disorder regime. Thus, having the condition $`\xi _1/L1`$ satisfied for each applied value of disorder, we can concentrate on the effects of disorder fluctuations by applying many configurations of random fields at fixed $`f`$ and $`L`$ values. We use up to $`L=768`$ and up to $`10^3`$ disorder configurations (i.e., up to $`36\times 10^6`$ spins).

Since we are interested in the properties of the triggering noise, which are related to distribution of disorder and the actual dynamics of the system, we avoid any “short-cuts” which can speed the algorithm . We apply the natural slow algorithm, which consists of the following steps: Random fields are stored at each site of an $`L\times L`$ lattice; The system is searched for the minimum local field $`h_{min}=\mathrm{min}(h_x^{loc})`$, where $`h_x^{loc}_x^{}J_{x,x^{}}S_x^{}+h_x+Hh_{ir,x}+H`$, and then the driving field is set to exactly $`h_{min}+\eta `$ (we use $`\eta =10^{10}`$), which thus triggers an avalanche at that site; A list is made consisting of the sites which may flip at the next time step (neighbor sites to flipped spins); The spins on the list are examined for flip and a new list is made; The process is continued until no more spins can flip, then next minimum local field is found. It should be stressed that all spins in a single list (spin shell) are updated in parallel, i.e., in a single time step, in analogy to parallel update in cellular automata models. In this way we have well defined time scale of avalanche evolution (internal time scale) as the number of steps that the updating procedure goes before the avalanche stops. On the other hand, the slow (external) time scale is set by the number of driving field updates.

It was shown in Ref.  that for low disorder $`f<f^{}=0.62`$ the extended domain wall moves through the system and depinns at a critical field $`H_c(f)`$, however, above the critical disorder $`f^{}`$ the built-in domain wall remains pinned and many new domains of reversed spins are nucleated inside the system. In this region, corresponding to high disorder (or low tensile stress), a typical structure of clusters which occurs in the linear part of the hysteresis is shown in Fig. 1. The theoretical value $`f^{}`$ can, in principle, be related to a critical value of tensile stress $`\sigma ^{}f_{}^{\mu }{}_{}{}^{}`$, below which the domain structure changes, as also observed in the experiment in Ref. . More precise experimental data would be necessary in order to determine the exponent $`\mu `$. A rough estimate is that $`\mu `$ is given by the correlation length exponent at the transition $`f^{}`$, $`\mu \nu 2.3`$ .

## III Fractal properties of triggering noise and domain-wall trapping

Time series of the magnetization jumps corresponding to the individual avalanches are known as BH noise. In the high disorder region either a new domain is nucleated and grown, or already existing domain extended, a BH pulse is associated with motion of a domain wall from one position to a new one. The area between two consecutive domain wall positions corresponds to the size of BH avalanche (measured as the area covered by a single BH pulse). In principle, a domain growth process is not “linear” but fractal, leading to the dynamic exponent $`z=1.23`$, and fractal dimension of avalanche $`D=1.88`$ . As usual, the dynamic exponent $`z`$ and the fractal dimension $`D`$ are defined via the scaling of characteristic duration and size of avalanches with the change of length scale, $`t_LL^z`$, and $`s_LL^D`$, respectively. Scaling properties are studied in detail (see set of scaling exponents in Ref. ). In particular, in the critical region above $`f^{}`$ the distribution of size of avalanches $`D(s)s^{\tau _s}`$ and duration $`P(t)t^{\tau _t}`$ scale with the exponents $`\tau _s=1.30`$ and $`\tau _t=1.47`$, respectively. The avalanche distributions integrated over hysteresis loop for $`f>f^{}`$ can be scaled according to the scaling form

$$P(t,f,L)=t^{\tau _t}𝒫(t(\delta f)^{z\nu },tL^z),$$
(2)
where $`\delta f(f/f^{}1)`$ and $`(\delta f)^\nu `$ measures the correlation length due to the critical point at $`f^{}`$. For instance, when the condition $`L(\delta f)^\nu `$ is satisfied, the avalanche cut-offs can be scaled with disorder strength (in the entire region of power-law behavior), leading to $`f^{}=0.62`$ and $`\nu =2.3`$ within statistical error bars (see  for details). It should be stressed that in this region of disorder the spatial extension of an avalanche is determined not only by strength of pinning but also by the blocking by previous clusters (see Fig. 1). The blocking effects dynamically alter the strength of pinning and thus influence the scaling properties. This illustrates a complex interplay between the disorder and dynamics, in contrast to, for instance, the equilibrium clusters in the random-field Ising model. In order to get an insight into these complex phenomena, we study next the distributions of wall trapping times and distances between the points with weakest pinning, from which an avalanche is released. We find that both of these distributions exhibit a scaling behavior.

The distribution of distances between initial points of the successive avalanches is shown in Fig. 2, where we distinguish between distances measured in the parallel ($`x_{}`$) and transverse() direction relative to the direction of the wall. In fact, due to anisotropy of domain wall motion in this case, these distributions show different scaling exponents $`\tau _{}=1.04\pm 0.03`$ and $`\tau _{}=0.60\pm 0.03`$ according to

$$G(x_{},x_{})x_{}^\tau _{}𝒢(x_{}/x_{}^{\zeta _G}),$$
(3)
where $`\zeta _G=\tau _{}/\tau _{}`$. Notice that for small distances correlations in both directions are almost equivalent. However, a crossover to a distinct power-law behavior for $`G(x_{})`$ occurs at a distance ($`x30`$ in Fig. 2) which is presumably related to the characteristic size of avalanche at given disorder. The open boundary in the direction opposite to wall leads to the cut-off at large $`x_{}`$. It should be stressed that spatial correlations are sensitive to distribution of disorder. Therefore, here we used a large number of samples in order to minimize scatter of the data due to disorder fluctuations and to clearly distinguish between correlations in parallel and perpendicular directions.

In the inset to Fig. 2 we show the Fourier spectrum of the threshold field fluctuations $`\mathrm{\Delta }H(t)`$ measured on the external time scale. It exhibits two correlated regions corresponding to earlier and later times, respectively. Steeper slope in the inset to Fig. 2 vary from 0.6-1.03, and flatter slope from 0.25 to 0.36, depending on $`f`$. The long-range correlations exhibited in Fig. 1 show that the system evolves in such way that its next relaxation event depends on the history of the present state of the system. It is interesting to note that a similar time series of the fluctuation of magnetization $`\mathrm{\Delta }M(t)M(t+1)M(t)`$ represents the Barkhausen noise itself (looked at the external time scale). An example of the BH noise signal is shown in Fig. 3. Note that on the external time scale duration of each elementary signal is equal to one, whereas the height of each elementary signal represents corresponding avalanche size. Sizes of the successive avalanches are only weakly correlated in time. However, a time derivative of the signal, namely $`a(t)d[\mathrm{\Delta }M(t)]/dt`$, representing acceleration of a domain wall at each field update, shows certain correlation properties. In the inset to Fig. 3 we show the Fourier spectrum of the numerical derivative of the the signal, which is shown in main Fig. 3. Average slope of the curve is close to one.

It should be noted that the average number of avalanches occurring at fixed strength of disorder increases with system size as $`n_aL^2/sL^{[2D(2\tau _s)]}\varphi (L(\delta f)^\nu )`$, where $`\varphi (L(\delta f)^\nu )`$ is unknown scaling function. This implies that in the theoretical case of infinitely slow driving the average size of the field jump $`\mathrm{\Delta }H(t)`$ decreases with $`L`$ as $`\mathrm{\Delta }H(t)H_{sat}(L)/n_aL^{0.68}H_{sat}(L)/\varphi (L(\delta f)^\nu )`$, where $`H_{sat}`$ is the saturation field. In real experimets the size of field jumps are restricted by the driving frequency, however, number of avalanches detected varies with the size of pick-up coil. This problem requires more detailed theoretical and experimental investigation .

Further understanding of the role of disorder in the dynamics can be achieved by considering the time intervals that a domain wall resides at a given point in space (trapping time $`T_{trap}`$). We calculate the distribution of trapping times $`T_{trap}`$ of a domain wall at a given site, which is determined as the time interval since the domain wall is pinned at a site $`(x_{},x_{})`$ until it eventually moves away from that site. In this way, $`T_{trap}`$ measures the time interval between two successive activities at that site. Notice that $`T_{trap}`$ in this case is somewhat different from so called first-return time in 1-dimensional interface . Here we have 2-dimensional system with many interacting interfaces. Another similar example is trapping of grains in rice-pile models . A reasonable time scale to measure $`T_{trap}`$ is the number of field updates (external time scale), since in the zero-temperature dynamics during an avalanche a spin at a given site is either fixed or flips only once.

The trapping-time distribution is shown in the inset to Fig. 4 for various values of disorder $`f>f^{}`$. Between the lower cut-off $`T_0`$ (below which all trapping times are equally probable), and an upper cut-off due to lattice size $`L`$, the trapping time distribution shows a short region with correlations. In particular, for large trapping times we can determine the slope as $`\tau _{trap}=2.3\pm 0.1`$. One can understand that a metastable configuration of domain walls corresponds to the system residing in a “pocket” of the fractal free energy landscape. A domain wall may move from a site $`(x_{},x_{})`$ as a part of an avalanche which started in the neighborhood of that point, corresponding to a local reorganization of the landscape near a shallow minimum (for a given value of the external field). When the system resides in a deeper minimum, however, it waits for larger driving fields (i.e., trapping time increases) or for a more global reconstruction of the landscape, which occurs with smaller probability compared to the one discussed above. This leads to the large slope of the distribution $`P(T_{trap})T_{trap}^{\tau _{trap}}`$ at large $`T_{trap}`$. It is interesting to note that a similar slope for large trapping times was found in the case of rice-pile model , where trapping of grains are considered, and it was argued that the slope is related to the roughness exponent of the rice-pile surface (see detailed analysis in Ref. ). Analogies between an interface motion and rice-pile model have been established in the literature . Considering the avalanche exponents the analogy also applies for the BH noise with an extended domain wall in the limit of low disorder . Notice that, in contrast to ricepile model, the range of power-law behavior in BH noise depends not only on the system size, but also on the strength of disorder, which restricts spatial extension of avalanches to $`ss_{max}L^2`$ and their durations to $`TT_{max}`$. Therefore, in this region of disorder we have rather small range of the power-law behavior. However, the scale-free behavior of the distribution of trapping times can be demonstrated via finite size scaling analysis when different system sizes are used. We find that the following scaling form

$$P(T_{trap},L)=L^\alpha 𝒫(T_{trap}L^{z_T}),$$
(4)
applies with the exponents which are weakly dependent on disorder $`\alpha =0.4\pm 0.05`$ and $`z_T=1.66\pm 0.05`$ (see Fig. 4). By increasing the system size $`L`$ the trapping times increase, leading to larger cut-offs of the distribution. We notice that the lower cut-offs $`T_0`$ also increase, leaving a limited correlation range. However, the cut-offs scale nicely with $`L`$, as shown in Fig. 4. In the inset to Fig. 4 we show how the trapping times distribution varies with disorder. By increasing disorder the cut-off $`T_0`$ moves towards larger values and the correlation region shrinks, corresponding to a lesser correlations in the dynamics, which also manifests in the decrease of the characteristic avalanche size. In the limit of infinitely strong disorder the dynamics becomes completely random, consisting of individual spin flips which align along a local random field.

It is interesting to note that unlike the avalanche distributions in Eq. (2), the distribution of trapping times in Eq. (4) as well as $`G(x)`$ (3) do not show explicit disorder dependence apart from a weak dependence in the exponents. In fact, the disorder effects in this case are included in the gradient of driving force $`\mathrm{\Delta }F`$ (i.e., the increments of the external field $`\mathrm{\Delta }H`$, which are adjusted to the minimum local fields). For example, for the distribution of trapping times we have in general

$$P(T_{trap},L,\mathrm{\Delta }F)=L^\alpha 𝒫(T_{trap}(\mathrm{\Delta }F)^{z_T/\lambda _F},T_{trap}L^{z_T}),$$
(5)
where the increase in the driving force $`\mathrm{\Delta }FH(t+T_{trap})H(t)`$ contains all the contributions due to interaction and random pinning which occurred during the time interval $`(t,t+T_{trap})`$, i.e., $`\mathrm{\Delta }F=_t^{t+T_{trap}}h_{ir,x}`$, which steam from the different sites in the system. According to the above results in Fig. 2, sum of these contributions has no characteristic scale. Therefore we may conclude that $`(\mathrm{\Delta }F)^{1/\lambda _F}\mathrm{}`$. Hence the first argument in the right hand side of Eq. (5) can be neglected compared to $`T_{trap}/L^{z_{trap}}`$, leading to the scaling form (4), which is in agreement with the scaling plot in Fig. 4.

Another interesting observation regards the dynamic exponent $`z_T`$, which scales the tail of the trapping time distribution with the system size $`L`$. It may be related to the scaling of the length of the optimal path  in strongly disordered medium, $`z_T=D_{OP}`$. In fact, the optimal path between two points can be constructed from the most persistent sections of the domain walls in the dynamics of BH avalanches. The length of the optimal path scales with the linear distance between the end points with the exponent $`D_{OP}=\tau _{OP}=1.66`$ .

## IV Transport equation with noise correlations

The long-range noise correlations are shown to be relevant for the scaling properties of the interface motion. In the literature the noise correlations in the interface depinning are viewed as originating from another (external) dynamic processes . In the case of Barkhausen avalanches we see that the noise correlations appear as an intrinsic property of the dynamics when the system is driven infinitely slowly. The active section of domain wall can be represented by a surface $`h(x,\tau )`$ which is pinned by quenched defects. The transport at the site $`(x,\tau )`$ is described by the equation

$$dh/d\tau =\nu _{}_{}^2h+\nu _{}_{}^2h+\eta (x,t,h).$$
(6)
Here we distinguish total elapsed time $`\tau `$, which is defined as the sum of evolution times of all individual avalanches, on one hand, and $`t`$, which represents time scale of field updates, on the other. The quenched noise $`\eta (x,t,h)`$ is generated by varying the driving field in the steps which are adjusted to the minimum local strength of pinning, as discussed above. We expect that a dominate $`h`$-dependence of the noise $`\eta `$ can be expressed as $`\eta p(x)h`$, where $`p(x)`$ represents the (anisotropic) local velocity of domain wall per field rate. In addition, time variation of the noise are related to the updates of the external field, and thus $`\eta `$ varies on the external time scale only. Therefore, we can write

$$\eta (x,t,h)p(x)(a_0+\mathrm{\Delta }H(t))(h),$$
(7)
where $`a_0`$ is a constant which depends on initial configuration and the hysteresis loop properties. The local interface velocity per field rate, $`p(x)`$, is a random variable which is determined by the spatial distribution and strength of pinning associated with a given value of the driving field $`H`$, i.e., history of the system. Therefore $`p(x)`$ is governed by a probability distribution, which can be deduced from the properties of the probability that an avalanche starts at distance $`x`$ from the preceding avalanche, i.e.,

$$<p(x)p(x^{})>=\gamma G(xx^{}).$$
(8)
It is assumed that the same type of correlations apply for the successive activities during the evolution of an avalanche, while the external field is constant. In the previous section we found that the probability $`G(x)`$ exhibits long-range correlations in the perpendicular and parallel direction as $`G(x_{},x_{})x_{}^1x_{}^{0.60}`$. We also notice that $`\mathrm{\Delta }H(t)t^\theta `$, where $`\theta `$ is finite in the case of infinitely slow driving discussed above, whereas, $`\theta =0`$ in the case of finite driving rate in small steps , where $`\mathrm{\Delta }H(t)=const`$. Notice that Eq. (6) together with Eqs. (7) and (8) leads to an effective nonlinear term of the form $`\gamma G(x_{},x_{})(a_0+t^\theta )^2(h)^2`$, which differs from the Kardar-Parisi-Zhang nonlinearity  by the power-law correlations in the coefficient. As discussed above, the origin of these correlations lies in the the dynamically varying pinning and blocking effects when a multidomain structure is slowly driven through the hysteresis loop.

The relevance of the noise correlations for the universality class of the interface depinning has been discussed in the literature using dynamic renormalization group (RG) . Another example where correlations of random noise play an important role is represented by the scaling properties of river networks . In that case the correlations of the type $`Gx_{}^{2+\delta }x_{}^{1\zeta }`$ were assumed, where $`\delta `$ is an expansion parameter and the anisotropy exponent $`\zeta `$ has to be determined self-consistently at a fixed point of the RG . In the result, the scaling exponents do depend on the range of correlations $`\delta `$. Similarly, we can expect that the long-range correlations of the form given in Eqs. (7) and (8) can be related to the universal scaling exponents of the BH noise. The exponents $`z`$, $`\zeta _G`$, which are defined above, and the roughness exponent $`\chi `$, which governs behavior of the dynamic variable $`h`$ with the change of scale, can in principle be obtained by the dynamic RG applied to Eq. (6) with noise properties in Eqs. (7-8). Then the avalanche exponents $`\tau _t`$ and $`\tau _s`$ can be deduced using scaling relations in the critical region. The RG analysis of the transport equation (6) with noise properties specified by Eqs. (7-8) requires additional work and is left out of the scope of this paper. Here we discuss only the scaling relations between the exponents at RG fixed point, e.g., the dynamic exponent $`z`$, and the avalanche exponents.

It should be stressed that such scaling relations are not universal and that they depend on the nature of the dynamic process (see for instance  for the case of river networks). We argue that the scaling relation derived below are valid for BH avalanches in the high-disorder (or multidomain) region. The rational behind these scaling relations is found in the directed nature of the avalanche propagation. The evolution of an avalanche (see Fig. 1) can be visualized as growth of a directed percolation (DP) cluster projected  to $`1+1`$ dimensions, with extra dimension representing time axis. The fractal dimension of the equivalent DP cluster measured with respect to the time axis is then $`D_{}D/z`$, where $`D`$ and $`z`$ are fractal dimension of BH avalanche and dynamic exponent, respectively, as defined above. Then $`D_{}=1+\zeta _{DP}\delta _{DP}`$, where $`\zeta _{DP}=z/2`$ is the anisotropy exponent measured with respect to time axis, and $`\delta _{DP}`$ is the survival probability exponent of the equivalent DP clusters. Notice also that for the directed dynamic processes $`\tau _t=D_{}`$ and thus $`\tau _s=21/\tau _t`$ , which completes the set of avalanche exponents. In addition, the roughness exponent $`\chi `$ can be related to the trapping time distribution as $`\tau _{trap}=2+\chi `$, as noticed in Refs. . For instance, using the numerical values for $`D`$ and $`z`$ from Ref. , the following values of the avalanche exponents are predicted by the above scaling relations: $`\tau _t=1.52`$, $`\tau _s=1.34`$, compared with $`1.47`$ and $`1.30`$, and $`\chi =0.23`$, obtained by direct numerical simulations.

## V Conclusions and discussions

We have demonstrated numerically the existence of long-range correlations in triggering noise which is intrinsic to the domain wall dynamics in slowly driven disordered ferromagnets. Although the distribution of disorder is initially uncorrelated, the spatio-temporal correlations develop in time. The appearance of these correlations can be attributed to the applied global driving in which always next weakest pinning force in the system is selected, on one hand, and to a finite extension of avalanches occurring between two consecutive field updates, on the other. Hence, attempts to reduce the domain wall dynamics to the Glauber spin flips in the presence of local random fields appears to be inadequate for the range of disorder where the cooperative avalanche dynamics occurs. We suggest an alternative transport equation which incorporates the observed noise correlations. The long-range correlations of triggering noise can be related to the universal scaling behavior of Barkhausen avalanches and to the transport properties via the fractal distribution of the domain wall trapping times. We believe that an analysis of the transport equation by the dynamic RG will contribute to understanding of the universality classes of Barkhausen avalanches by the infinitely slow driving. By varying the driving conditions, however, these correlations are changed. This might be the origin of different scaling behavior of BH avalanches at finite driving rates, as observed in experiments . Our results also suggest that technical procedures which speed the numeric algorithms in driven disordered systems should be taken with care. In particular, an algorithm which alters the properties of triggering noise may lead to a different scaling behavior, which is unrelated with the original problem.

###### Acknowledgements.

This work was supported by the Ministry of Science and Technology of the Republic of Slovenia. I am grateful to Deepak Dhar for many helpful comments and suggestions.


format: text

text:
# Spectral properties and pseudogaps in a model with d-wave pairing symmetry

\[

## Abstract

A model with d-wave pairing symmetry is studied by employing a non-perturbative sum rule approach. At low temperature the magnitude of a normal state pseudogap shows strong $`\stackrel{}{k}`$ or angle dependence well fitted by $`\mathrm{cos}2\varphi `$ form. With increasing temperature, the pseudogap closes at some critical angle $`\varphi _c`$ and beyond this angle a single quasiparticle-like peak appears. The resulting Fermi surface is strongly temperature dependent. Both in the spectral function and the density of states, the pseudogap disappears in a manner that the spectral weight fills in the pseudogap instead of closing it with increasing temperature. All these features are qualitatively consistent with ARPES for underdoped cuprates.

\]

The nature of an excitation gap in high-temperature superconductors has been one of the puzzling issues in the community of condensed matter physicists. Through several years of extensive experimental work, general consensus regarding the superconducting gap symmetry seems to be reached that the superconducting gap has mainly d-wave character with possibility of a small mixture of other angular momentum states, in contrast to conventional BCS superconductors with an isotropic s-wave gap. Recent discovery of a normal state pseudogap in underdoped cuprates has shed another side of the anomalous behaviors in the copper oxide superconductors. For these materials the low frequency spectral weight begins to be strongly suppressed below some characteristic temperature $`T^{}`$ higher than $`T_c`$. This behavior has been observed through various experimental probes such as photoemission, specific heat, tunneling, NMR, and optical conductivity. In particular recent angle resolved photoemission spectroscopy (ARPES) and tunneling experiments indicate that the pseudogap phenomenon is closely related to pairing fluctuations. These measurements clearly exhibit that the normal state pseudogap has the same angular dependence and magnitude as the superconducting gap and that often the only difference between the spectra in the pseudogap state and the superconducting state is in their linewidths. Typically $`T^{}`$ is much higher than $`T_c`$ and their doping dependence is qualitatively different. While $`T_c`$ decreases with underdoping, $`T^{}`$ increases in contrast. This feature suggests that $`T^{}`$ does not follow $`T_c`$ characterized by long-range phase coherence, but instead some kind of a mean-field critical temperature $`T_{MF}`$. In spite of several possible scenarios such as the spinon pair formation without the Bose-Einstein condensation of holons, strong superconducting phase fluctuations and a magnetic scenario near the antiferromagnetic instability and so on, at present there is no consensus in the origin of the pseudogap.

For the past several years extensive theoretical effort has been also made by several groups to understand this anomalous pseudogap behavior in the context of short range effective (attractive) interaction between electrons. This may be divided into two different classes of approach. In the first class, quantum Monte Carlo (QMC) simulations were made for the attractive Hubbard model. Although this Hamiltonian is not a realistic model for understanding the complex physical behaviors in the underdoped cuprates, it is believed to capture the important ingredient of the paring fluctuations in those materials. In the absence of a small parameter this numerical method has played an important role in understanding of the model, in spite of some uncertainties due to finite size effect and numerical analytical continuation. In this approach the one-particle spectral function as well as the density of states clearly show the precursor of the superconducting gap in the normal state. In the second class, various assumptions and approximations are used. This class includes the effect of vortex phase fluctuations on the single-particle properties and the ‘paring approximation’ theory, and T-matrix and self-consistent T-matrix approximations for the two-dimensional attractive Hubbard model. More recently the self-consistent T-matrix approach was also applied to a model with $`d_{x^2y^2}`$ pairing. The present approach is a variant of the T-matrix approximation applied to a model with d-wave pairing symmetry.

In a previous study, a non-perturbative sum rule approach was developed for the attractive Hubbard model by extending previous work on the repulsive Hubbard model. It is found that in two dimensions, the mean-field transition temperature is replaced by a crossover temperature where the correlation length starts to grow exponentially. At sufficiently low temperature, a Kosterlitz-Thouless O(2) transition should occur, but it is not reproduced by that approach since it is in the $`O(n=\mathrm{})`$ universality class. Nevertheless, the agreement with Monte Carlo calculations is quantitative for both one and two-particle correlation functions over the whole range of parameters accessible by Monte Carlo calculations where it is found that in two dimensions, a pseudogap appears in the single-particle spectral weight as well as in the density of states. Recent dimensional crossover study by Preosti et al. shows that the pseudogap effect is basically absent in three dimensions. In weak to intermediate coupling, the appearance of a pseudogap is traced back to the growing critical pairing fluctuations in the low-temperature renormalized classical regime of the low-dimensional system. With increasing temperature, the spectral weight fills in the pseudogap instead of closing it. Furthermore, the pseudogap appears earlier in the density of states than in the spectral function. It was also noted that the qualitative features found in this study should apply to the d-wave case. In this paper we study in detail spectral properties and pseudogaps in a model Hamiltonian with d-wave pairing symmetry which is more appropriate for high-temperature superconductors.

We consider on a two dimensional square lattice a simple model Hamiltonian which has a superconducting ground state with d-wave symmetry

$`H`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{k},\sigma }{}}\epsilon _\stackrel{}{k}c_{\stackrel{}{k},\sigma }^+c_{\stackrel{}{k},\sigma }`$ (1)
$`+`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{}{k},\stackrel{}{k}^{},\stackrel{}{q}}{}}V_{\stackrel{}{k},\stackrel{}{k}^{}}c_{\stackrel{}{k},}^+c_{\stackrel{}{k}+\stackrel{}{q},}^+c_{\stackrel{}{k}^{}+\stackrel{}{q},}c_{\stackrel{}{k}^{},},`$ (2)
where $`\epsilon _\stackrel{}{k}=2t(\mathrm{cos}k_x+\mathrm{cos}k_y)`$ and $`N`$ is the number of lattice sites. We choose $`V_{\stackrel{}{k},\stackrel{}{k}^{}}`$ as a d-wave separable potential given as $`V\mathrm{\Pi }(\stackrel{}{k})\mathrm{\Pi }(\stackrel{}{k}^{})`$ where $`\mathrm{\Pi }(\stackrel{}{k})=\mathrm{cos}k_x\mathrm{cos}k_y`$, or as $`1.5264\mathrm{cos}2\varphi =1.5264\frac{k_x^2k_y^2}{k_x^2+k_y^2}`$ where $`\varphi =\mathrm{arctan}\frac{k_y}{k_x}`$. In the real space notation the interaction term may be written as $`H_\text{I}=V_i\mathrm{\Delta }_i^+\mathrm{\Delta }_i`$ where $`\mathrm{\Delta }_i=_\delta g(\delta )c_{i,}c_{i+\delta ,}`$. $`g(\delta )`$ is the Fourier transform of $`\mathrm{\Pi }(\stackrel{}{k})`$ given as $`g(\delta )=\frac{1}{N}_\stackrel{}{k}e^{i\stackrel{}{k}\stackrel{}{\delta }}\mathrm{\Pi }(\stackrel{}{k})`$. For example, for the following pair structure,

$`g(\delta )=\{\begin{array}{ccc}1/2\hfill & \text{if }\delta =(\pm 1,0),\hfill & \\ 1/2\hfill & \text{if }\delta =(0,\pm 1),\hfill & \\ 0\hfill & \text{if otherwise},\hfill & \end{array}`$

it leads to $`\mathrm{\Pi }(\stackrel{}{k})=\mathrm{cos}k_x\mathrm{cos}k_y`$. Obviously $`\mathrm{\Pi }(\stackrel{}{k})=1.5264\mathrm{cos}2\varphi `$ requires in $`g(\delta )`$ further neighbors beyond the four nearest neighbors. Here a numerical factor of 1.5264 is introduced to normalize $`g(\delta )`$ with $`_\delta g^2(\delta )=1`$. In terms of $`\mathrm{\Delta }_i`$ the d-wave pairing susceptibility is defined as

$`\chi _{pp}(i,\tau )=T_\tau \mathrm{\Delta }_i(\tau )\mathrm{\Delta }_0^+,`$ (3)
where $`T_\tau `$ is the imaginary time ordering operator.

It is well-known that for many strongly correlated electronic models there is no obvious small parameter with which a systematic perturbative approximation can be made. In such a situation, a non-perturbative approach may be a good alternative for at least qualitative understanding of the physics in those strongly correlated model Hamiltonians. In this paper we apply a non-perturbative sum rule approach to the model Hamiltonian (2), which is still not a fully controlled approximation in nature. Because of an extended nature of $`\mathrm{\Delta }_i`$ for a d-wave pair and the increasing difficulty of the derivation in the real space representation, here we use the analogy from our previous studies for the repulsive and attractive Hubbard models. These studies show the important many-body modification with respect to the standard RPA and T-matrix approaches. The modification comes in two different places. First, the paring susceptibility is calculated through vertex function $`U_{pp}`$ instead of bare interaction strength $`U`$, which is constant in our approximation. Second, the pairing fluctuations are coupled to electrons, leading to the self-energy with $`UU_{pp}`$ form instead of $`U^2`$. In particular, the latter structure is in agreement with the fact that there is no Migdal theorem for this problem, contrary to the case of electron-phonon interactions. And the renormalized interaction constant $`U_{pp}`$ is determined by the exact sum rule for the pairing susceptibility. By using the sum rule and the renormalized constant $`U_{pp}`$, we determine in effect the Ginzburg-Landau parameters due to mode-mode coupling. Thus this approach is similar to the one-loop renormalization group approach within the Gaussian approximation and as a result the Mermin-Wager theorem is formally satisfied in two dimensions. This sum rule approach was systematically compared with the QMC simulations for the repulsive and attractive Hubbard models and the agreement was in a quantitative level both in one- and two-particle functions. We believe that this generic feature of the many-body modification should also carry over to a Hamiltonian with d-wave pairing symmetry. In this modification, the pairing susceptibility and the self-energy can be written in Fourier space in terms of renormalized interaction strength $`V_{pp}`$

$`\chi _{pp}(q)`$ $`=`$ $`{\displaystyle \frac{\chi _{pp}^0(q)}{1+V_{pp}\chi _{pp}^0(q)}},`$ (4)
$`\mathrm{\Sigma }(k)`$ $`=`$ $`VV_{pp}\mathrm{\Pi }^2(\stackrel{}{k}){\displaystyle \frac{T}{N}}{\displaystyle \underset{q}{}}\chi _{pp}(q)G^0(qk),`$ (5)
where the irreducible susceptibility is defined as

$`\chi _{pp}^0(q)={\displaystyle \frac{T}{N}}{\displaystyle \underset{k}{}}\mathrm{\Pi }^2(\stackrel{}{k})G^0(qk)G^0(k).`$ (6)
The above expressions are reduced to the standard T-matrix approximation when $`V_{pp}`$ is replaced by bare $`V`$. We determine this renormalized constant by employing the exact sum rule for the d-wave pairing susceptibility

$`{\displaystyle \frac{T}{N}}{\displaystyle \underset{q}{}}\chi _{pp}(q)e^{i\nu _m0^{}}`$ (7)
$`=`$ $`{\displaystyle \underset{\delta ,\delta ^{}}{}}g(\delta )g(\delta ^{})c_{\delta ^{},}^+c_{\delta ,}c_{}^+c_{}.`$ (8)
Note that the convergence factor is necessary in the sum rule, because away from half-filling $`\chi _{pp}(i\nu _m)`$ decays like $`1/\nu _m`$ at large frequencies. In the previous studies for the repulsive and attractive Hubbard models, the right-hand side of the sum rule evaluated in the SDW and BCS mean-field ground states, respectively, was found to be in excellent agreement with the QMC values in the intermediate to strong coupling regimes. In this paper the right-hand side of Eq. (8) is also evaluated in the d-wave BCS mean-field ground state

$`{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{\stackrel{}{k},\stackrel{}{p},\stackrel{}{k}^{},\stackrel{}{p}^{}}{}}\mathrm{\Pi }(\stackrel{}{k})\mathrm{\Pi }(\stackrel{}{p})c_{\stackrel{}{k},}^+c_{\stackrel{}{p},}c_{\stackrel{}{k}^{},}^+c_{\stackrel{}{p}^{},}\delta _{\stackrel{}{k}+\stackrel{}{k}^{},\stackrel{}{p}+\stackrel{}{p}^{}}`$ (9)
$`=`$ $`[{\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{}{k}}{}}\mathrm{\Pi }^2(\stackrel{}{k})v^2(\stackrel{}{k})][{\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{}{k}}{}}v^2(\stackrel{}{k})]`$ (10)
$`+`$ $`[{\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{}{k}}{}}\mathrm{\Pi }(\stackrel{}{k})u(\stackrel{}{k})v(\stackrel{}{k})]^2,`$ (11)
where

$`u_\stackrel{}{k}^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1+{\displaystyle \frac{\epsilon _\stackrel{}{k}\mu }{E_\stackrel{}{k}}}),`$ (12)
$`v_\stackrel{}{k}^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1{\displaystyle \frac{\epsilon _\stackrel{}{k}\mu }{E_\stackrel{}{k}}}),`$ (13)
$`E_\stackrel{}{k}`$ $`=`$ $`\sqrt{(\epsilon _\stackrel{}{k}\mu )^2+\mathrm{\Delta }^2(\stackrel{}{k})}.`$ (14)
Here $`\mathrm{\Delta }(\stackrel{}{k})`$ is the d-wave BCS mean-field gap. The chemical potential $`\mu `$ and the gap $`\mathrm{\Delta }(\stackrel{}{k})`$ are determined self-consistently through the number and gap equations for given $`V`$, $`T`$ and $`n`$.

Before starting we comment on some differences associated with $`\mathrm{\Pi }(\stackrel{}{k})=\mathrm{cos}k_x\mathrm{cos}k_y`$ and with $`\mathrm{\Pi }(\stackrel{}{k})=1.5264\mathrm{cos}2\varphi `$ structure. The interaction Hamiltonian with the first structure depends not only on the angle of $`\stackrel{}{k}`$ but also on its magnitude. Thus the pairing interaction is always strongest at $`\stackrel{}{k},\stackrel{}{k}^{}=(\pm \pi ,0)`$ or $`(0,\pm \pi )`$, although the noninteracting Fermi surface can be far from these points. As a result two features occur for a particle density far away from half-filling (Detailed calculations were performed but not shown in this paper). First, the peaks associated with the precursor of the Bogoliubov quasiparticles occur asymmetrically with respect to the Fermi energy. Second, the locus of $`\stackrel{}{k}`$ points satisfying $`\omega \epsilon _\stackrel{}{k}+\mu Re\mathrm{\Sigma }(\stackrel{}{k},\omega )=0`$ at $`\omega =0`$ can be substantially different from the noninteracting Fermi surface, thus strongly violating the Luttinger’s theorem. For $`\mathrm{\Pi }(\stackrel{}{k})=1.5264\mathrm{cos}2\varphi `$ which depends only on the angle, however, the above features disappear and the locus of $`\stackrel{}{k}`$ points satisfying the above equation is almost identical with the noninteracting Fermi surface. Near half-filling the differences are negligible. Throughout the calculations, lattice spacing, $`\mathrm{}`$, and $`k_B`$ are set to be unity, and all energies are measured in unit of $`t`$. We used a discrete lattice as large as $`128\times 128`$ in momentum space and performed the calculations by means of fast Fourier transforms (FFT). Equations (5) and (6) are computed in terms of Matsubara frequencies and the analytic continuation from Matsubara to real frequencies are made via Pade approximants. In order to detect any spurious features associated with numerical analytical continuation, we performed real frequency calculations in parallel. Except for some spiky features in the real frequency formulation coming from the Lorentzian approximation of the non-interacting Green’s function, the two results are almost identical.

We begin by presenting the spectral functions along $`(0,0)(\pi ,0)`$ direction for $`V=4t`$, $`n=0.5`$ and $`T=0.15t`$ in Fig. 1. Throughout the paper the density is fixed at $`n=0.5`$. The results for other densities are similar. For this parameter, $`V_{pp}`$ satisfying the sum rule is found to be $`1.87t`$ significantly different from the bare value ($`V=4t`$). This shows the importance of the mode-mode coupling effect already in the intermediate coupling regime. Let us label wave vectors by $`(m\pi /8,0)`$. Below the Fermi wave vector ($`m=5`$), the main peak stays below the Fermi energy and the secondary peak grows in strength as $`\stackrel{}{k}`$ approaches $`\stackrel{}{k}_F`$. At the Fermi wave vector two peaks reminiscent of the Bogoliubov quasiparticles appear almost symmetrically with respect to the Fermi energy. Since electrons are still in the normal state, this is the precursor of the superconducting gap, namely, a normal state pseudogap. Above the Fermi wave vector, the main peak stays above the Fermi energy and the secondary peak becomes stronger as $`\stackrel{}{k}`$ approaches $`\stackrel{}{k}_F`$. This result should be contrasted with that by Engelbrecht et al. In the self-consistent T-matrix approximation with $`\mathrm{\Pi }(\stackrel{}{k})=\mathrm{cos}k_x\mathrm{cos}k_y`$ form, these authors argued that along $`(0,0)(\pi ,0)`$ direction the dominant peak of $`A(\stackrel{}{k},\omega )`$ never crosses the Fermi energy and bounces back towards the negative frequency. Their overall finding, however, is qualitatively different from our results. In the present calculations, the dominant peak of $`A(\stackrel{}{k},\omega )`$ eventually passes through the Fermi energy for $`\stackrel{}{k}`$ far above the Fermi wave vector. At the noninteracting Fermi surface our spectral weight is strongly suppressed at the Fermi energy to become a local minimum, while the local minimum of $`A(\stackrel{}{k},\omega )`$ in their study is located significantly away from the Fermi energy. Presumably this is due to the features associated with $`\mathrm{\Pi }(\stackrel{}{k})=\mathrm{cos}k_x\mathrm{cos}k_y`$ structure mentioned in the previous paragraph and also due to the self-consistent approximation that does not take the vertex function and the Green’s function at the same level of approximation.

In Fig. 2 the spectral function is shown along the Fermi wave vectors at low temperature ($`T=0.15t`$). The angle $`\varphi `$ is defined as $`\mathrm{arctan}(k_y/k_x)`$ along the noninteracting Fermi surface. At $`\varphi =0^o`$ the magnitude of the pseudogap is largest and as $`\varphi `$ increases it progressively decreases, leading to strong momentum or angle dependence in the size of the pseudogap. Since along the diagonal directions the pairing interaction vanishes, the pseudogap completely closes and a sharp quasiparticle peak appears at $`\varphi =45^o`$. In Fig. 3 the magnitude of the pseudogap (circles) is plotted as a function of angle $`\varphi `$ for $`T=0.15t`$ along with the ground state gap symmetry $`\mathrm{\Delta }(\varphi =0^o)\mathrm{cos}2\varphi `$ (dashed curve). The angle dependence of the normal state pseudogap is well fitted by d-wave symmetry, consistent with ARPES experiments for underdoped cuprates.

At higher temperatures a drastic change is found in the spectral function. For $`T=0.225t`$ the pseudogap closes well below $`45^o`$, as shown in Fig. 4. The local minimum at the Fermi energy disappears at $`\varphi =31.0^o`$ ($`26.6^o`$ is just on the crossover). Beyond this angle, the quasiparticle-like peak appears and thereby the Fermi surface is partially restored. We define this angle as critical angle $`\varphi _c`$. In Fig. 5 $`\varphi _c`$ is plotted for different temperatures (open circles). Below $`T=0.175t`$ the pseudogap is found everywhere along the Fermi wave vectors except at $`\varphi =45^o`$, thus the Fermi surface is destroyed everywhere except along the diagonal directions. With increasing temperature (up to $`T=0.29t`$), however, $`\varphi _c`$ becomes smaller than $`45^o`$. Thus, the pseudogap region shrinks and at the same time the Fermi surface grows from $`\varphi =45^o`$ up to the critical angle. Above $`T=0.29t`$ the whole Fermi surface is completely restored in spite of some broadening of the spectral function due to interaction as well as finite temperature. This feature is qualitatively consistent with ARPES for underdoped cuprates. We can theoretically calculate the temperature dependence of $`\varphi _c`$. By using the Ornstein-Zernike form of the pairing correlation function and taking the classical fluctuations ($`iq_n=0`$), the scattering rate at the Fermi energy is found to be proportional to $`\mathrm{\Pi }^2(\stackrel{}{k})\xi /\xi _T`$, where $`\xi `$ and $`\xi _T`$ are pairing correlation length and thermal de Broglie wave length defined as $`\xi _T=v_F(\stackrel{}{k})/T`$, respectively. For the pseudogap to disappear, the scattering rates should be much smaller than unity $`\mathrm{\Pi }^2(\stackrel{}{k})\xi /\xi _T1`$, allowing us to define a critical angle $`\xi =v_F(\varphi _c)/T/\mathrm{\Pi }^2(\varphi _c)`$. In Fig. 5 this critical angle is also shown as stars. For the best fit near $`0^o`$ angle, a numerical factor of 1.06 multiplies $`\xi `$. Although a small deviation is found near $`45^o`$, the overall magnitude and shape are in reasonable agreement throughout the whole angle. For a d-wave model, the condition for the appearance of a pseudogap in a given momentum depends not only on the anisotropy of the Fermi velocity (which is the only relevant condition in the attractive Hubbard model) but also more importantly on an angle dependent form factor $`\mathrm{\Pi }^2(\stackrel{}{k})`$.

In Fig. 6 we show the spectral function at the Fermi wave vector and the density of states for different temperatures. As the temperature is increased, the spectral weight starts to fill in the pseudogap and at $`T=0.3t`$ the precursor of the superconducting gap completely disappears as shown in Fig. 6(a). At this temperature, however, the pseudogap still persists in the density of states and at higher temperature ($`T=0.4t`$) it finally disappears as shown in Fig. 6(b). Except close to half-filling, the pseudogap appears at higher temperature in the density of states than in the spectral function. Compared with the s-wave case, the density of states is suppressed linearly near the Fermi energy, a reminiscence of $`N(\omega )\omega `$ in the superconducting state. Both in the spectral function and the density of states, the pseudogap disappears in a manner that the spectral weight fills in the pseudogap instead of closing it with increasing temperature. This feature is also consistent with ARPES for underdoped cuprates. This may suggest that in our approach phase fluctuations (spin-wave type) rather than amplitude fluctuations are mainly responsible for the pseudogap formation, although the present approach includes both. Like in the s-wave case, a pseudogap also appears in the density of states when the characteristic pairing frequency $`\nu _c`$ is equal to or smaller than temperature. This corroborates the origin of the d-wave pseudogap, namely, growing d-wave paring fluctuations in the low temperature renormalized classical regime of the low dimensional system. The calculated $`T^{}`$ follows the same trend as the mean-field critical temperature $`T_{MF}`$. (To be more precise, $`T^{}`$ is approximately half of $`T_{MF}`$ for most of the densities.) As noted in Ref., near a point with high order parameter symmetry (half-filling in that paper), the transition temperature $`T_c`$ decreases while the pseudogap temperature increases along with $`T_{MF}`$. As a result it leads to a large pseudogap regime, consistent with the phase diagram in the underdoped side of cuprates.

There are several advantages in our formulation. First, the pairing fluctuation sum rule Eq. (8) is exactly satisfied (by construction). Through this sum rule, the Mermin-Wagner theorem is formally fulfilled and the strength of pairing fluctuations is properly constrained within the Gaussian approximation. This latter feature is crucial in our formulation, because an approximate treatment of pairing fluctuations without constraining the strength can easily overestimate or underestimate the magnitude of fluctuations particularly in low dimensional systems. The Mermin-Wagner theorem is also satisfied in some other approaches such as the self-consistent T-matrix (FLEX) and the “pairing” approximation schemes. Second, there is an exact relation between one-particle (self-energy, Green’s function) and two-particle (interaction term) functions:

$`\underset{\tau 0^{}}{lim}{\displaystyle \frac{T}{N}}{\displaystyle \underset{\stackrel{}{k},i\omega _n}{}}\mathrm{\Sigma }(\stackrel{}{k},i\omega _n)G(\stackrel{}{k},i\omega _n)e^{i\omega _n\tau }=V\mathrm{\Delta }_i^+\mathrm{\Delta }_i,`$ (15)
where here the self-energy includes the Hartree-Fock term. When the interacting Green’s function $`G(\stackrel{}{k},i\omega _n)`$ is replaced by the noninteracting one $`G^0(\stackrel{}{k},i\omega _n)`$ in Eq. (15), it is exactly satisfied. With $`G(\stackrel{}{k},i\omega _n)`$, the difference between the left- and right-hand sides is less than $`6\%`$ for all temperatures studied. In this paper the Kosterlitz-Thouless phase transition and its fluctuation effect near $`T_{KT}=T_c`$ have been neglected, since the present formulation is inadequate to describe the topological nature of the phase in two dimensions and vortex anti-vortex binding-unbinding physics. Since for similar parameters for the attractive Hubbard model the QMC results indicate that $`T_{KT}0.05t`$ three times smaller than the lowest temperature in our calculations, we do not expect in the present results any significant influence from vortex phase fluctuations. Our approach is valid for weak to intermediate coupling and for temperature not too deep in the pseudogap regime. At very low temperature, even the Ginzburg-Landau functional form itself may change, for instance, a possible crossover of the dynamical critical exponent from z=2 to other value, and eventually vortex phase fluctuations may come into play.

In summary, we have studied spectral properties and pseudogaps in a model with d-wave pairing symmetry by using a non-perturbative sum rule approach. We applied to this model our previous experience of many-body theory in the repulsive and attractive Hubbard models. The magnitude of the normal state pseudogap shows strong angle dependence well fitted by $`\mathrm{cos}2\varphi `$ form at low temperature. With increasing temperature the pseudogap closes at some critical angle $`\varphi _c`$ and beyond this angle a single quasiparticle-like peak appears. The resulting Fermi surface is strongly temperature dependent. With increasing temperature, the pseudogap region shrinks and at the same time the Fermi arc grows from $`\varphi =45^o`$ to $`\varphi _c`$. Both in the spectral function and density of states the pseudogap disappears in a manner that the spectral weight fills in the pseudogap instead of closing it with increasing temperature. All these features are qualitatively consistent with ARPES for underdoped cuprates. The pseudogap is caused by growing d-wave critical pairing fluctuations in the low-temperature classical renormalized regime of the low-dimensional system, as in the repulsive and attractive Hubbard models. We argue that although the real critical behaviors and critical exponents are governed by the vortex phase fluctuations close to the $`T_{KT}`$, the initial growth of pairing fluctuations can be driven by spin-wave phase fluctuations, leading to the normal state pseudogap formation. The behavior of the spin-wave type phase fluctuations belonging to the $`O(n=\mathrm{})`$ universality class can be qualitatively different in two dimensions (particularly near a point with high order parameter symmetry) from that of the $`O(2)`$ Kosterlitz-Thouless vortex phase fluctuations.

The author would like to thank A. M. Tremblay for numerous help and discussions throughout the work. This work was supported by a grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada and the Fonds pour la formation de Chercheurs et d’Aide à la Recherche (FCAR) of the Québec government.