format: text

text:
# 1 Plot of the potentials 𝑉₂(𝑥) of (22) as a function of 𝑥 and ℰ when 𝑒^{2𝜆}=-1/(2ℰ) and 𝜈=0. The oscillator potential 𝑥²/2-1 is a member of this family arising for ℰ=-1/2.

ON FIRST-ORDER SCALING INTERTWINING IN QUANTUM MECHANICS

\[Rev. Mex. Fís. 46 S2, (Nov. 2000) 153-156\]

DAVID J. FERNÁNDEZ C.

Departamento de Física, CINVESTAV-IPN, Apdo Postal 14-740, 07000 México D.F., Mexico

HARET C. ROSU

Instituto de Física, Universidad de Guanajuato, Apdo Postal E-143, 37150 León, Gto, Mexico

Abstract. We generalize the standard first-order intertwining relationship of supersymmetric quantum mechanics in order to include simultaneous scaling transformations in both the original Hamiltonian and the intertwining operator. It is argued that in this way one can generate potentials with more interesting spectra than those obtained by means of the standard first-order intertwining technique and, as an outcome, a simple engineering procedure is presented. The harmonic oscillator potential is used in order to illustrate the previous statements. Moreover, a matrix representation of the scaled intertwining relationship is sketched up allowing for higher-dimensional generalizations in the case of separable potentials.
Resumen. Generalizamos la relación de entrelazamiento estandar de primer orden de la mecánica cuántica supersimétrica para incluir de manera simultánea transformaciones de escalamiento tanto en el Hamiltoniano original como en el operador de entrelazamiento. Se argumenta que en esta forma uno puede generar potenciales con espectros más interesantes que aquellos obtenidos por medio de la técnica de entrelazamiento estandar de primer orden y, como un resultado, un procedimiento sencillo de ingenieria cuántica es presentado. El potencial de oscilador armónico es usado para ilustrar las afirmaciones anteriores. Mas aún, una representación matricial de la relación de entrelazamiento escalado es bosquejado que permite la generalización a más de una dimensión en el caso de los potenciales separables.

Factorizations of second order linear one-dimensional (1D) differential operators are common tools in Witten’s supersymmetric quantum mechanics (SUSYQM) \[?\], which may be considered as a form of Darboux transformations \[?\]. They imply either particular solutions of Riccati equations known as superpotentials or the general Riccati solution (RS), the latter case being first used in physics by Mielnik for the quantum harmonic oscillator \[?\]. In recent years, it became clear that SUSYQM is not only a form of Darboux transformations but also a particular case in the more general framework provided by the technique of intertwining operators \[?\], of extensive use in the mathematical literature. Two (Hamiltonian) operators are said to be intertwined by an operator $`A^+`$ if the following relationship is fulfilled

$$H_1A^+=A^+H.$$
(1)
In SUSYQM $`A^+`$ is a first order differential (factorization) operator of the form

$$A^+=\frac{1}{\sqrt{2}}\left(\frac{d}{dx}+\alpha (x)\right),$$
(2)
leading to the standard Riccati equation for $`\alpha (x)`$ associated to the given initial potential

$$\alpha ^{}(x)+\alpha ^2(x)=2(V(x)),$$
(3)
where the prime denotes derivative with respect to $`x`$. The potential $`V_1(x)`$ corresponding to $`H_1`$ is determined according to:

$$V_1(x)=V(x)\alpha ^{}(x),$$
(4)
whenever one is able to find a solution of (3) for given $`V(x)`$ and $``$. The so-called factorization energy $`𝐑`$ plays a crucial role in generating new solvable potentials from a given one. This becomes clear when substituting $`\alpha (x)=[\mathrm{ln}\psi (x)]^{}`$ in (3), which leads to

$$H\psi (x)=\frac{1}{2}\psi ^{\prime \prime }(x)+V(x)\psi (x)=\psi (x).$$
(5)
Although similar to the standard eigenvalue equation for $`H`$, notice that $`\psi (x)`$ in (5) is not necessarily normalizable, but it should not have zeros in order to avoid supplementary singularities of $`V_1(x)`$ with respect to those of $`V(x)`$. It is well known that for $`>E_0`$, where $`E_0`$ is the ground state energy of $`H`$, $`\psi (x)`$ will always have zeros. However, if $`E_0`$ it is possible to make $`\psi (x)`$ to have no zeros by adjusting the ratio of the two constants in the general solution of (5), resulting in a physically meaningful $`V_1(x)`$ as explained for example by Sukumar \[?\]. The spectrum of $`V_1(x)`$ consists of the sequence of levels $`E_n`$ of the initial Hamiltonian plus a new ‘ground state’ energy level at $``$. The scheme becomes complete after realizing that $`H`$ and $`H_1`$ are factorized in terms of $`A`$ and $`A^+`$ as follows:

$$H=AA^++,H_1=A^+A+.$$
(6)
If now one writes $`H_1=Hf(x)`$, the intertwining relationship (1) can be put in the commutator form as follows

$$[H,A^+]=f(x)A^+,$$
(7)
and moreover, it is easy to see that Eq. (4) requires

$$f(x)=\alpha ^{}(x),$$
(8)
that is, in the standard SUSYQM the function $`f(x)`$ is always the derivative of the RS $`\alpha (x)`$, and it will be called the Darboux potential difference (DPD) between the Darboux pair of almost isospectral potentials $`V(x)`$ and $`V_1(x)`$, respectively.

We now look for a more general first-order intertwining relationship and perform two simultaneous modifications: one related to the Hamiltonian in which we introduce a “perturbative” parameter $`ϵ`$, $`H_1=(1+ϵ)Hf(x)`$, and the other related to the intertwining operator (2) which is multiplied (say, to the right) by a unitary operator $`Ue^{i\frac{\lambda }{2}(xp+px)}`$ ($`U^+=U^1`$) depending on a ‘scaling’ real parameter $`\lambda `$:

$$A_\lambda ^+=\frac{1}{\sqrt{2}}\left(D+\alpha (x)\right)U,$$
(9)
where $`D\frac{d}{dx}`$. Multiplying Eq. (1) by $`U^+`$ to the right one gets

$$[(1+ϵ)Hf(x)](D+\alpha (x))=(D+\alpha (x))UHU^+,$$
(10)
and making equal on both sides of (10) the coefficients of the various corresponding powers of the derivative operator $`D`$ leads to

$`D^3:`$ $`e^{2\lambda }=1+ϵ,`$ (11)
$`D^2:`$ $`e^{2\lambda }\alpha (x)=(1+ϵ)\alpha (x),`$ (12)
$`D^1:`$ $`(1+ϵ)\alpha ^{}(x)+(1+ϵ)V(x)f(x)=V(e^\lambda x),`$ (13)
$`D^0:`$ $`{\displaystyle \frac{(1+ϵ)}{2}}\alpha ^{\prime \prime }(x)+[(1+ϵ)V(x)f(x)]\alpha (x)=\alpha (x)V(e^\lambda x)V^{}(e^\lambda x).`$ (14)
One can see that Eqs. (11) and (12) are the same, giving the relationship between the parameters $`ϵ`$ and $`\lambda `$. After substituting $`e^{2\lambda }=1+ϵ`$, Eq. (13) reads

$$f(x)=e^{2\lambda }\alpha ^{}(x)+e^{2\lambda }V(x)V(e^\lambda x),$$
(15)
which is a more general form of Eq. (8). We can see that the original potential at two different points interfere in the connection between the DPD $`f(x)`$ and the RS $`\alpha (x)`$. This makes a key difference between the ‘scaling’ intertwining and the standard SUSYQM one, even for the cases when the potentials are homogeneous functions of degree $`d`$ for which (15) simplifies further to the more local form $`f(x)=s^2\alpha ^{}(x)+(s^2s^d)V(x)`$, where $`s=e^\lambda `$. When $`d=2`$ one gets $`f(x)=e^{2\lambda }\alpha ^{}(x)`$ which is the closest one can reach to the standard SUSYQM. Substituting Eqs. (11) and (15) in Eq. (14) and integrating it we get the Riccati equation

$$e^{2\lambda }\alpha ^{}(x)+e^{2\lambda }\alpha ^2(x)=2(V(e^\lambda x)),$$
(16)
where, as it will be clear below, it is convenient to choose the integration constant equal to the factorization energy $``$ of (3). Now, if one makes the change of variable $`y=e^\lambda x`$ and denotes $`\stackrel{~}{\alpha }(y)=e^\lambda \alpha (e^\lambda y)`$ one easily gets

$$\frac{d\stackrel{~}{\alpha }(y)}{dy}+\stackrel{~}{\alpha }^2(y)=2(V(y)).$$
(17)
One can see that, up to a rescaling of the coordinate, this is the Riccati equation (3) with a factorization energy $``$. Therefore, the same RS used in the standard intertwining technique (Eqs. (1) to (8)) can be used as well in the generalized scaling intertwining technique of Eqs.(1,9-17) in order to generate new solvable potentials with known spectra. The eigenfunctions of $`H_2e^{2\lambda }H_1=He^{2\lambda }f(x)`$ are proportional to the action of $`A_\lambda ^+`$ on the eigenfunctions of $`H`$ due to the fact that the two Schrödinger Hamiltonians $`H`$ and $`H_2`$ are also intertwined, although now in a slightly generalized way:

$$H_2A_\lambda ^+=e^{2\lambda }A_\lambda ^+H.$$
(18)
The potential of the scaled intertwined Hamiltonian $`H_2`$ takes the form:

$$V_2(x)=V(x)e^{2\lambda }f(x)=e^{2\lambda }V(e^\lambda x)\alpha ^{}(x).$$
(19)
The corresponding eigenvalues are $`\{e^{2\lambda },e^{2\lambda }E_n\}`$, where $`e^{2\lambda }`$ is the ground state energy level of $`H_2`$ associated to the eigenfunction $`\psi _{}(x)\mathrm{exp}(_0^x\alpha (z)𝑑z)`$. The factorizations of $`H`$ and $`H_2`$ in terms of $`A_\lambda `$ and $`A_\lambda ^+`$ become:

$$H=e^{2\lambda }A_\lambda A_\lambda ^++,H_2=A_\lambda ^+A_\lambda +e^{2\lambda }.$$
(20)
The generation of new exactly solvable potentials $`V_2(x)`$ by means of the scaling intertwining technique depends on the expertness to solve the Riccati equation (17). In particular, if either we find the general solution for just an isolated value of the factorization energy $``$ or fixing it if we have the solution for $``$ in a real interval, the family of potentials $`V_2(x)`$ derived by means of this technique would have the spectrum $`\{e^{2\lambda },e^{2\lambda }E_n\}`$, and by varying $`\lambda `$ we would be scaling the basic spectrum $`\{,E_n\}`$ at $`\lambda =0`$ generated by means of the standard intertwining technique when creating a new level at $``$ starting from the initial potential $`V(x)`$. However, the really interesting procedure is as follows. Suppose we would solve (17) for $``$ belonging to a real interval $`I=[_1,_2]`$, where $`_1`$ and $`_2`$ have the same sign. By choosing now $`=e^{2\lambda }_0`$ with $`_0I`$ fixed and varying $`\lambda `$ such that always $`I`$, we would obtain potentials $`V_2(x)`$ having spectra of the type $`\{_0,(\frac{_0}{})E_n\}`$, i.e., changed so that the excited state levels would be scaled while the ground state energy level of $`H_2`$ would be unaffected. This interesting effect is different of the one achieved by using the standard intertwining technique, where after solving the Riccati equation (3) for $``$ belonging to a real interval, when varying $``$ in its corresponding domain we would change the ground state energy level but maintaining unaffected the excited state levels. Thus, these two kinds of intertwining effects, alone or combined, provide us with considerable freedom for designing solvable potentials with suitable spectra for modeling purposes.

As an example, let us consider the harmonic oscillator potential $`V(x)=x^2/2`$. The solution to the Riccati equation (17) for an arbitrary $`<1/2`$ is given by (see, e.g., Sukumar and Junker and Roy ):

$$\stackrel{~}{\alpha }(y)=y+\frac{d}{dy}\{\mathrm{ln}\left[{}_{1}{}^{}F_{1}^{}(\frac{1+2}{4},\frac{1}{2};y^2)+2\nu \frac{\mathrm{\Gamma }(\frac{32}{4})}{\mathrm{\Gamma }(\frac{12}{4})}y{}_{1}{}^{}F_{1}^{}(\frac{3+2}{4},\frac{3}{2};y^2)\right]\},$$
(21)
where, in order to avoid singularities we have that $`|\nu |<1`$. Thus, the 3-parametric family of potentials (the parameters are $`\lambda ,,\nu `$),

$$V_2(x)=e^{4\lambda }\frac{x^2}{2}e^\lambda \stackrel{~}{\alpha }^{}(e^\lambda x),$$
(22)
with $`\stackrel{~}{\alpha }(y)`$ given by (21), has spectrum $`\{e^{2\lambda },e^{2\lambda }(n+1/2),n=0,1,\mathrm{}\}`$. By varying $`\lambda `$ and maintaining $``$ fixed, we will move into members of the family (22) attainable after applying the scale transformation onto the potentials generated by means of the standard intertwining technique, which have spectrum $`\{,E_n=n+1/2,n=0,1,\mathrm{}\}`$. On the other hand, if we take $`=e^{2\lambda }_0I_1=(\mathrm{},0)`$, where $`_0I_1`$ is fixed, when varying the scaling parameter $`\lambda `$ we would be changing $``$ so that one will get a two-parameter family of potentials with spectrum $`\{_0,(\frac{_0}{})(n+1/2),n=0,1,\mathrm{}\}`$, i.e., the excited state levels will be scaled by the factor $`_0/`$, but the ground state energy level will remain fixed. A similar treatment can be implemented for $`I_2=(0,1/2)`$. A plot of the potentials $`V_2(x)`$ as a function of $`x`$ and $`I_1`$, with $`_0=1/2`$ and $`\nu =0`$, is shown in Figure 1. Notice that this selection of $`_0`$ and $`\nu `$ ensures that, up to a displacement of the energy origin, the oscillator potential is included in this family (for $`=1/2`$, i.e., $`\lambda =0`$), and for other values of $``$ $`V_2(x)`$ is symmetric with respect to $`x=0`$. For $`(1/2,0)`$ the potentials $`V_2(x)`$ have a double well, and the corresponding spacing for the excited state levels is expanded by the factor $`1/(2)`$ ($`1`$ is the original spacing of the oscillator potential). The ground state energy level in this case is fixed at $`1/2`$. On the other hand, for $`(\mathrm{},1/2)`$, the potentials $`V_2(x)`$ present just a peaked single well centered at $`x=0`$, which can be seen as a deformation of the well of the oscillator potential $`x^2/21`$ needed to maintain fixed the ground state energy level at $`1/2`$. The excited state levels are ‘squeezed’ by the factor $`1/(2)`$.

In interesting works, Spiridonov \[?\] has used a rather similar type of scaling intertwining in a different notation (the $`q`$ one) for emphasizing the connection between the $`q`$-deformed calculus \[?\] and Shabat’s infinite chain of reflectionless potentials \[?\]. In such a treatment the so-called self-similar potentials have been introduced and characterized, and at the same time the corresponding spectrum has been determined by means of the algebraic technique itself. On the other hand, in this paper the orientation has been different: to take initially a potential with known spectrum in order to generate a new potential with known and different (in general) spectrum. Since in the limit $`\lambda 0`$ we recover the results derived by means of the standard first order intertwining, we have at hand an appropriate and interesting generalization of such a technique. For any other $`\lambda `$ the new spectrum will not have common elements with the initial one, and moreover, the above mentioned non-local influence of the initial potential on the DPD (see Eq. (15)) will be at work. This is one of the interesting novelties put forth by means of scaling/deformation within the SUSYQM intertwining procedure.

The generalization of the scaled intertwining formalism to more dimensions can be easily performed for the cases when one deals with separable potentials by means of the matrix representation \[?\]. Thus, for separable 2D potentials one introduces $`𝒜^{}=\left(\begin{array}{cc}A& 0\\ 0& B\end{array}\right)`$ and $`𝒜^+=\left(\begin{array}{cc}A^+& 0\\ 0& B^+\end{array}\right)`$. The symbols $`A`$ and $`B`$ denote the factorization operators for the first and second coordinate axis, respectively. The 2$`\times `$2 matrix intertwining relationship corresponding to (1) reads

$$_1𝒜^+=𝒜^+,$$
(23)
where

$$=\left(\begin{array}{cc}AA^++_1& 0\\ 0& BB^++_2\end{array}\right).$$
(24)
Corresponding to the function $`f(x)`$ and the unitary operator $`U`$ one introduces the following 2$`\times `$2 diagonal matrices

$$=\left(\begin{array}{cc}f_1(x)& 0\\ 0& f_2(y)\end{array}\right)=\left(\begin{array}{cc}\alpha _1^{}(x)& 0\\ 0& \alpha _2^{}(y)\end{array}\right)$$
(25)
and

$$𝒰=\left(\begin{array}{cc}e^{i\frac{\lambda _1}{2}(xp_x+p_xx)}& 0\\ 0& e^{i\frac{\lambda _2}{2}(yp_y+p_yy)}\end{array}\right),$$
(26)
respectively. In (26) the two scaling parameters may or may not be equal. By means of these matrices, all the 1D formulas of this paper have 2$`\times `$2 matrix counterparts, though one can immediately see that in order to maintain the diagonality of matrices one needs separable potentials. As an example, the matrix form of equation (7) reads $`[,𝒜^+]=𝒜^+`$. It is also quite easy to generalize this matrix representation to any number of even dimensions.

Finally, we draw attention to the fact that the scaling intertwining formalism could be applied to any type of Hamiltonians displaying scaling properties, not necessarily quantum ones. On the other hand, many interesting 2D cases are not separable. We recall as an example the diamagnetic Kepler problem, i.e., the Hamiltonian problem of the hydrogen atom in a magnetic field \[?\], that in atomic units ($`\mathrm{}=e=m_e=1`$) and for zero angular momentum reads $`H=\frac{1}{2}𝐩^2+\frac{1}{8}\gamma ^2(x^2+y^2)\frac{1}{r}`$, where $`\gamma =B`$/(2.35$`\times `$$`10^5`$T) gives the magnetic field strength. In scaled coordinates and momenta, $`\stackrel{~}{𝐫}=\gamma ^{2/3}𝐫`$ and $`\stackrel{~}{𝐩}=\gamma ^{1/3}𝐩`$, respectively, this Hamiltonian can be written $`\stackrel{~}{H}=\gamma ^{2/3}H`$ and does not depend on the energy and field strength separately, but only on the scaled energy $`\stackrel{~}{E}=\gamma ^{2/3}E`$. This case, which is also of much laboratory and astrophysical interest, despite its scaling features, might be suitable for a more general scaling intertwining approach with coupled Riccati equations and non-diagonal matrices. Since here our purposes are merely illustrative we refer the reader to our recent preprint \[?\] for more complicated cases.

Acknowledgements

This work was supported in part by the CONACyT Projects 26329-E and 458100-5-25844E. HCR wishes to thank for the kind hospitality at the Departamento de Física, CINVESTAV.

References


text:
# Transverse energy distributions and 𝐽/𝜓 production in Pb+Pb collisions

## Abstract

We have analyzed the latest NA50 data on transverse energy distributions and $`J/\psi `$ suppression in Pb+Pb collisions. The transverse energy distribution was analysed in the geometric model of AA collisions. In the geometric model, fluctuations in the number of NN collisions at fixed impact parameter are taken into account. Analysis suggests that in Pb+Pb collisions, individual NN collisions produces less $`<E_T>`$, than in other AA collisions. The nucleons are more transparent in Pb+Pb collisions. The transverse energy dependence of the $`J/\psi `$ suppression was obtained following the model of Blaizot et al, where charmonium suppression is assumed to be 100% effective above a threshold density. With fluctuations in number of NN collisions taken into account, good fit to the data is obtained, with a single parameter, the threshold density.

PACS numbers: 25.75.-q, 25.75.Dw

In relativistic heavy ion collisions $`J/\psi `$ suppression has been recognized as an important tool to identify the possible phase transition to quark-gluon plasma. Because of the large mass of the charm quarks, $`c\overline{c}`$ pairs are produced on a short time scale. Their tight binding also make them immune to final state interactions. Their evolution probes the state of matter in the early stage of the collisions. Matsui and Satz  predicted that in presence of quark-gluon plasma (QGP), binding of $`c\overline{c}`$ pairs into $`J/\psi `$ meson will be hindered, leading to the so called $`J/\psi `$ suppression in heavy ion collisions . Over the years several groups have measured the $`J/\psi `$ yield in heavy ion collisions (for a review of the data and the interpretations see ). In brief, experimental data do show suppression. However this could be attributed to the conventional nuclear absorption, also present in pA collisions.

The latest data obtained by the NA50 collaboration  on J/$`\psi `$ production in Pb+Pb collisions at 158 A GeV is the first indication of anomalous mechanism of charmonium suppression, which goes beyond the conventional suppression in nuclear environment. The ratio of J/$`\psi `$ yield to that of Drell-Yan pairs decreases faster with $`E_T`$ in the most central collisions than in the less central ones. It has been suggested that the resulting pattern can be understood in a deconfinement scenario in terms of successive melting of charmonium bound states .

In a recent paper Blaizot et al  showed that the data can be understood as an effect of transverse energy fluctuations in central heavy ion collisions. Introducing a factor $`\epsilon =E_T/E_T(b)`$ and assuming that the suppression is 100% above a threshold density (a parameter in the model) and smearing the threshold density (at the expense of another parameter) best fit to the data was obtained. Capella et al  analysed the data in the comover approach. There also, the comover density has to be modified by the factor $`\epsilon `$. Introduction of this adhoc factor $`\epsilon `$ can be justified in a model based on excited nucleons represented by strings .

At a fixed impact parameter, the transverse energy as well as the number of NN collisions fluctuate. The Fluctuations in the number of NN collisions were not taken into account in the calculations of Blaizot et al  or in the calculations of Capella et al . In the present paper, we present a calculation following the model of Blaizot et al  which includes these fluctuations. As will be shown below, if the fluctuations in number of NN collisions are taken into account very good description of the NA50 data can be obtained without smearing the threshold density. The smearing effect is generated by the fluctuations.

Geometric model has been quite successful in explaining the transverse energy production as well as multiplicity distributions in AA collisions . In this model, it is assumed that at impact parameter $`𝐛`$, the number $`n`$ of NN collisions is Poisson distributed with average $`<n_b>`$. In the Glauber approximation $`<n_b>`$ is written as,,

$$<n_b>=\sigma _{NN}d^2sT_A(𝐬)T_B(𝐬𝐛)$$
(1)
where $`\sigma _{NN}`$ is the inelastic NN cross-section, assumed to be 32 mb. In this model all the nuclear information is contained in the nuclear thickness function is $`T_{A,B}(𝐬)=𝑑z\rho _{A,B}(𝐬,z)`$. In the present calculation, we have used the following parametric form for $`\rho _A(r)`$ ,

$$\rho _A(r)=\frac{\rho _0}{1+exp(\frac{rr_0}{a})}$$
(2)
with a=0.53 fm, $`r_0=1.1A^{1/3}`$. The central density is obtained from $`\rho _A(r)d^3r=A`$.

In the geometric model, the probability to obtain $`E_T`$ at impact parameter $`𝐛`$ in $`n`$\- number of NN collisions is written as,

$$P_n(b,E_T)=\frac{e^{<n_b>}<n_b>^n}{\mathrm{\Gamma }(n+1)}Q^{\{n\}}(E_T)$$
(3)
where $`Q^{\{n\}}(E_T)`$ is the n-fold convolution of $`E_T`$ distribution resulting from elementary NN collisions,

$$Q^{\{n\}}(E_T)=𝑑E_T^1\mathrm{}𝑑E_T^ng(E_T^1)\mathrm{}g(E_T^n)\times \delta (E_TE_T^1\mathrm{}E_T^n)$$
(4)
In eq.4 $`g(E_T)`$ is the normalized $`E_T`$ distribution for NN collisions. Most of the $`E_T`$ distributions in NN collisions can be well approximated by the gamma distribution,

$$g(E_T)=\frac{\alpha ^\beta }{\mathrm{\Gamma }(\beta )}e^{\alpha E_T}E_T^{\beta 1}$$
(5)
with the parameters $`\alpha `$ and $`\beta `$. For the gamma distribution, the average and the variance are,

$`<E_T>_{NN}=`$ $`\beta /\alpha `$ (7)
$`<E_T^2>/<E_T>^21=`$ $`1/\beta `$ (8)
Gamma distribution has an elegant convolution property which greatly facilitate computation. n-fold convolution of a gamma distribution is again a gamma distribution, with parameters $`\alpha ^{}=\alpha `$ and $`\beta ^{}=n\beta `$. Thus,

$$Q^{\{n\}}(E_T)=\frac{\alpha ^{n\beta }}{\mathrm{\Gamma }(n\beta )}e^{\alpha E_T}E_T^{n\beta 1}$$
(9)
The final transverse energy distribution is then obtained from eq.3 by summing it over n (from 1 to $`\mathrm{}`$) and averaging over the impact parameter b.

Geometric model has been quite successful in explaining the $`E_T`$ distributions in heavy ion collisions with $`\alpha `$ 2 and $`\beta `$ 2 . We have fitted the $`E_T`$ distribution in Pb+Pb collisions varying the parameter $`\alpha `$ and $`\beta `$. NA50 collaboration did not correct the $`E_T`$ spectra for the efficiency of target identification algorithm, which is lower than unity for $`E_T`$ lower than 60 GeV. To obtain the parameters $`\alpha `$ and $`\beta `$ we have fitted the purely inclusive part of the $`E_T`$ spectra ($`E_T>`$ 60 GeV). Very good fit to the data is obtained with $`\alpha =3.46\pm .19`$ and $`\beta =0.379\pm 0.021`$. The fit is shown in fig.1. We note that best fitted values of $`\alpha `$ and $`\beta `$ indicate that average $`E_T`$ produced in individual NN collisions is order of magnitude smaller than in other AA collisions . The variance is also a order of magnitude large. It seems that the nucleons are more transparent in Pb+Pb collisions than in (say) O+Au collisions, leading to less $`<E_T>`$ in individual NN collisions. It is evident that $`E_T`$ production mechanism in Pb+Pb collisions is different from other AA collisions. Transparency of nucleons in Pb+Pb collisions may be interpreted as an indication of QGP production in Pb+Pb collisions.

As mentioned in the beginning, we have followed the model of Blaizot et al  to analyze the $`E_T`$ dependence of charmonium suppression. Charmonium production cross-section at impact parameter $`𝐛`$ is written as,

$$d^2\sigma ^{J/\psi }/d^2b=\sigma ^{J/\psi }d^2sT_A^{eff}(𝐬)T_B^{eff}(𝐬𝐛)S(𝐛,𝐬)$$
(10)
where $`T_{A,B}^{eff}`$ is the effective nuclear thickness function,

$$T^{eff}(𝐬)=_{\mathrm{}}^{\mathrm{}}𝑑z\rho (𝐬,z)exp(\sigma _{abs}_z^{\mathrm{}}𝑑z^{}\rho (𝐬,z^{}))$$
(11)
with $`\sigma _{abs}`$ as the cross-section for $`J/\psi `$ absorption by nucleons. The exponential factor is the nuclear absorption survival probability, the probability for the $`c\overline{c}`$ pair to avoid nuclear absorption and form a $`J/\psi `$. $`S(𝐛,𝐬)`$ is the anomalous part of the suppression. Blaizot et al  assumed that $`J/\psi `$ suppression is 100% effective above a threshold density ($`n_c`$), a parameter in the model. Accordingly the anomalous suppression part was written as,

$$S(𝐛,𝐬)=\mathrm{\Theta }(n_cn_p(𝐛,𝐬))$$
(12)
where $`n_p`$ is the density of participant nucleons in impact parameter space,

$$n_p(𝐛,𝐬)=T_A(𝐬)[1e^{\sigma _{NN}T_B(𝐛𝐬)}]+[T_AT_B]$$
(13)
Recognizing that the endpoint behavior of charmonium suppressions are due to transverse energy fluctuations, Blaizot et al  modified the density of participant nucleons by a factor $`\epsilon =E_T/E_T(b)`$. This modification makes sense only when $`n_p`$ is assumed to be proportional to the energy density. Implicitly it was also assumed that the $`E_T`$ fluctuations are strongly correlated in different rapidity gaps. The assumption was essential as NA50 collaboration measured $`E_T`$ in the 1.1-2.3 pseudorapidity window while the $`J/\psi `$’s were measured in the rapidity window $`2.82<y<3.92`$ . In the geometric model, strong correlation between $`E_T`$ fluctuations in different rapidity window is readily obtained. At a fixed impact parameter, fluctuations in $`E_T`$ can be calculated as

$$\frac{<E_T^2>_{AA}<E_T>_{AA}^2}{<E_T>_{AA}^2}=\frac{<N^2><N>^2}{<N>^2}+\frac{1}{<N>}\frac{<E_T^2>_{NN}<E_T>_{NN}^2}{<E_T>_{NN}^2}$$
(14)
The $`E_T`$ fluctuations has two parts, (i) geometric in nature which remains same irrespective of rapidity window, (ii) which depend on fluctuations in the $`E_T`$ distributions in NN collisions. The 2nd part changes with rapidity window but its effect is less as it is weighted by the factor $`1/<N>`$. Fluctuations of $`E_T`$ in different rapidity windows are thus strongly correlated.

We calculate the $`J/\psi `$ production as a function of transverse energy, at an impact parameter $`𝐛`$ as,

$$d\sigma ^{J/\psi }/dE_T=\underset{n=1}{\overset{\mathrm{}}{}}P_n(b,E_T)P(\psi E_T,𝐛)$$
(15)
where $`P_n(b,E_T)`$ is the probability to obtain $`E_T`$ in $`n`$ NN collisions (eq.(3)) and $`P(\psi E_T,𝐛)`$ is the probability to produce a charmonium with transverse energy $`E_T`$. $`P(\psi E_T,𝐛)`$ is given by eq.10, with anomalous suppression part modified according to,

$$S(𝐛,𝐬)=\mathrm{\Theta }(n_c\frac{E_T}{n\beta /\alpha }n_p(𝐛,𝐬))$$
(16)
where we have replaced $`E_T(b)`$ by $`n\beta /\alpha `$ appropriate in the geometric model. This modification takes into account the fluctuations in number of NN collisions at fixed impact parameter $`𝐛`$.

The Drell-Yan production was calculated similarly, replacing $`P(\psi E_T,𝐛)`$ in eq.15 by the Drell-Yan production cross-section,

$$d^2\sigma ^{DY}/d^2b=\sigma ^{DY}d^2sT_A(𝐬)T_B(𝐬𝐛)$$
(17)
In fig.2, we have compared the theoretical charmonium production cross-section with NA50 experimental data. The normalization factor $`\sigma ^{J/\psi }/\sigma ^{DY}`$ was taken to be 53.5. The solid curve is obtained with $`\sigma _{abs}`$=6.4 mb, and $`n_c=3.8fm^2`$. Very good description of the data from 40 GeV onward is obtained. It may be noted that if the fluctuations in the NN collisions were neglected, equivalent description is obtained with threshold density $`n_c=3.75fm^2`$, with smearing of the $`\mathrm{\Theta }`$ function at the expense of another parameter. It is evident that in this model, the smearing is done by fluctuating NN collisions. Theoretical calculations predict more suppressions below 40 GeV, a feature evident in other models also. It is possible to fit the entire $`E_T`$ range, reducing the $`J/\psi `$-nucleon absorption cross-scetion. The dashed line in fig.2, corresponds to $`\sigma _{abs}`$=4 mb and $`n_c=3.42fm^2`$.

To summarize, we have analysed the transverse energy distribution in Pb+Pb collisions as well as the charmonium production data obtained by the NA50 collaboration. The transverse energy distribution was analysed in the geometric model. It was seen that in order to fit the experimental data, individual NN collisions are required to produce (on the average) less $`E_T`$, compared to other AA collisions. It seems that in Pb+Pb collision, the nucleons become transparent. The charmonium production data was analysed following the model of Blaizot et al , including the effect of fluctuations in number of NN collisions at fixed impact parameter. The experimental data from 40 GeV onwards could be very well fitted with $`\sigma _{abs}`$=6.4 mb and a threshold density of 3.8 $`fm^2`$. Neglecting the fluctuations in number of NN collisions, equivalent fit could only be obtained by smearing the $`\mathrm{\Theta }`$ function at the expense of an added parameter. Data in the entire $`E_T`$ range could be fitted reducing the $`J/\psi `$-nucleon absorption cross-section to $`\sigma _{abs}`$=4 mb. The threshold density is also reduced to 3.44 $`fm^2`$. Considering that nucleons become more transparent (as suggested by the $`E_T`$ data), such a reduction seems plausible. Melting of charmoniums above a threshold density as well as apparent transparency of nucleons (as evident from analysis of $`E_T`$ distribution) strongly suggests that in Pb+Pb collisions, QGP like environment is produced.


text:
# Quasilinearization Approach to Nonlinear Problems in Physics with Application to Nonlinear ODEs

## I Introduction

In a series of recent papers,  the possibility of applying a very powerful approximation technique called the quasilinearization method (QLM) to physical problems has been discussed. The QLM is designed to confront the nonlinear aspects of physical processes. The method, whose iterations are constructed to yield rapid convergence and often monotonicity, was originally introduced forty years ago by Bellman and Kalaba  as a generalization of the Newton-Raphson method  to solve individual or systems of nonlinear ordinary and partial differential equations. Modern developments and applications of the QLM to different fields are given in a monograph .

However, the QLM was never systematically studied or extensively applied in physics, although references to it can be found in well known monographs  dealing with the variable phase approach to potential scattering, as well as in a few scattered research papers . The reason for the sparse use of the QLM in Physics is that the convergence of the method has been proven only under rather restrictive conditions , which generally are not fulfilled in physical applications. Recently, though, it was shown  by one of the present authors (VBM) that a different proof of the convergence can be provided which we will generalize and elaborate here so that the applicability of the method is extended to incorporate realistic physical conditions of forces defined on infinite intervals with possible singularities at certain points.

In the first paper of the series , the quasilinearization approach was applied to the nonlinear Calogero equation in a variable phase approach to quantum mechanics and the results were compared with those of perturbation theory and the exact solutions. It was found analytically and by examples that the $`n`$-th approximation of the QLM exactly sums $`2^n1`$ terms of perturbation theory. In addition, a similar number of terms is reproduced approximately. The number of the exactly reproduced perturbation terms thus doubles with each subsequent QLM approximation, and reaches, for example, 127 terms in the 6-th QLM approximation, 8191 terms in the 12-th QLM approximation, and so on.

The computational approach in the work  was mostly analytical, and therefore one was able to compute only two to three QLM iterations, mainly for power potentials. Only in the case of the $`1/r^2`$ potential, could the calculation of QLM iterations be done analytically for any $`n`$.

The goal of the next work  was, by dropping the restriction of analytical computation, to calculate higher iterations as well as to extend the analysis to non-power potentials, in order to better assess the applicability of the method and of its numerical stability and the convergence pattern of the QLM iterations. It was shown that the first few iterations already provide very accurate and numerically stable answers for small and intermediate values of the coupling constant and that the number of iterations necessary to reach a given precision only moderately increases for larger values of the coupling. The method provided accurate and stable answers for any coupling strengths, including for super singular potentials for which each term of the perturbation theory diverges and the perturbation expansion does not exist even for a very small coupling.

The quasilinearization approach is applicable to a general nonlinear ordinary or partial $`n`$ -th order differential equation in $`N`$-dimensional space. In this paper, we consider the case of nonlinear ordinary differential equations in one variable which, unlike the nonlinear Calogero equation considered in references, contain not only quadratic nonlinear terms but various other forms of nonlinearity and not only a first, but also higher derivatives. Namely, we apply it to a panopoly of well-known and difficult nonlinear ordinary first, second and third order differential equations and show that again with just a small number of iterations one can obtained fast convergent and uniformly excellent and stable numerical results.

The paper is arranged as follows: in the second chapter we present the main features of the quasilinearization approach, while in the third chapter we consider, as a warm-up exercise, a simple first-order differential equation with a nonlinear $`n`$-th power term and compare its exact analytic solution with the perturbation theory and with the QLM iterations in order to demonstrate the main features of the quasilinearization approach. In the next four chapters, we apply our method to four well known nonlinear ordinary second and third order differential equations, namely to the Lane-Emden, Thomas-Fermi, Duffing, and Blasius equations, respectively. These equations have been extensively studied in the literature. For a sample of recent papers see refs.  and references therein.

The results, convergence patterns, numerical stability, advantages of the method and its possible future applications are discussed in the final chapter.

## II The quasilinearization method (QLM)

The aim of the quasilinearization method (QLM) of Bellman and Kalaba  based on the Newton-Raphson method  is to solve a nonlinear $`n`$-th order ordinary or partial differential equation in $`N`$ dimensions as a limit of a sequence of linear differential equations. This goal is easily understandable since there is no useful technique for obtaining the general solution of a nonlinear equation in terms of a finite set of particular solutions, in contrast to a linear equation which can often be solved analytically or numerically in a convenient fashion using superposition. In addition, the QL sequence should be constructed to assure quadratic convergence and, if possible, monotonicity.

As we have mentioned in the Introduction, we will follow here the derivation outlined in ref. , which is not based, unlike the derivations in refs. , on a smallness of the interval and on the boundness of the nonlinear term and its functional derivatives, the conditions which usually are not fulfilled in physical applications.

For simplicity, we limit our discussion to nonlinear ordinary differential equation in one variable on the interval $`[0,b],`$ which could be infinite:

$$Ł^{(n)}u(x)=f(u(x),u^{(1)}(x),\mathrm{}.u^{(n1)}(x),x),$$
(1)
with $`n`$ boundary conditions

$$g_k(u(0),u^{(1)}(0),\mathrm{}.,u^{(n1)}(0))=0,k=1,\mathrm{}l$$
(2)
and

$$g_k(u(b),u^{(1)}(b),\mathrm{}.,u^{(n1)}(b))=0,k=l+1,\mathrm{},n.$$
(3)
Here $`L^{(n)}`$ is linear $`n`$-th order ordinary differential operator and $`f`$ and $`g_1,g_2,\mathrm{}..,g_n`$ are nonlinear functions of $`u(x)`$ and its $`n1`$ derivatives $`u^{(s)}(x),s=1,\mathrm{}n1`$. The more general case of partial differential equations in N-dimensional space could be considered in exactly the same fashion by changing the definition of $`L^{(n)}`$ to be a linear $`n`$-th order differential operator in partial derivatives and $`x`$ to be an N-dimensional coordinate array.

The QLM prescription  determines the $`r+1`$-th iterative approximation $`u_{r+1}(x)`$ to the solution of Eq. (1) as a solution of the linear differential equation

$`L^{(n)}u_{r+1}(x)=f(u_r(x),u_r^{(1)}(x),\mathrm{}..,u_r^{(n1)}(x),x)`$ (4)
$`+{\displaystyle \underset{s=0}{\overset{n1}{}}}(u_{r+1}^{(s)}(x)u_r^{(s)}(x))f_{u^{(s)}}(u_r(x),u_r^{(1)}(x),\mathrm{}..,u_r^{(n1)}(x),x),`$ (5)
where $`u_r^{(0)}(x)=u_r(x)`$, with linearized two-point boundary conditions

$$\underset{s=0}{\overset{n1}{}}(u_{r+1}^{(s)}(0)u_r^{(s)}(0))g_{ku^{(s)}}(u_r(0),u_r^{(1)}(0),\mathrm{}..,u_r^{(n1)}(0),0)=0,k=1,\mathrm{}l$$
(6)
and

$$\underset{s=0}{\overset{n1}{}}(u_{r+1}^{(s)}(b)u_r^{(s)}(b))g_{ku^{(s)}}(u_r(b),u_r^{(1)}(b),\mathrm{}..,u_r^{(n1)}(b),b)=0,k=l+1,\mathrm{},n.$$
(7)
Here the functions $`f_{u^{(s)}}=f/u^{(s)}`$ and $`g_{ku^{(s)}}=g_k/u^{(s)},s=0,1,\mathrm{},n1`$ are functional derivatives of the functionals $`f(u(x),u^{(1)}(x),\mathrm{}.u^{(n1)}(x),x)`$ and $`g_k(u(x),u^{(1)}(x),\mathrm{}.u^{(n1)}(x),x)`$, respectively. For example, in case of a simple nonlinear boundary condition $`u^{}(b)u(b)=c`$ where c is a constant, one has $`g(r)g(u(r),u^{}(r),r)=u^{}(r)u(r)`$ so that $`g_u=u^{}(r)`$ and $`g_u^{}=u(r)`$. The linearized boundary condition 7 has a form $`(u_{r+1}(b)u_r(b))u_r^{}(b)+(u_{r+1}^{}(b)u_r^{}(b))u(b)=0`$ or $`(u_{r+1}(b)u_r(b))^{}=(u_r(b)u_r(b))^{}`$ so the nonlinear boundary condition for the initial guess $`u_0(b)u_0^{}(b)=c`$ will be propagated to the linear boundary condition for the next iterations..

The zeroth approximation $`u_0(x)`$ is chosen from mathematical or physical considerations.

To prove that the above procedure yields a quadratic and often monotonic convergence to the solution of Eq. 1 with the boundary conditions 2 and 3, we follow reference  and consider a differential equation for the difference $`\delta u_{r+1}(x)u_{r+1}(x)u_r(x)`$ between two subsequent iterations:

$`L^{(n)}\delta u_{r+1}(x)`$ $`=`$ $`[f(u_r(x),u_r^{(1)}(x),\mathrm{}..,u_r^{(n1)}(x),x)f(u_{r1}(x),u_{r1}^{(1)}(x),\mathrm{}..,u_{r1}^{(n1)}(x),x)]`$ (8)
$`+`$ $`{\displaystyle \underset{s=0}{\overset{n1}{}}}[\delta u_{r+1}^{(s)}(x)f_{u^{(s)}}(u_r(x),u_r^{(1)}(x),\mathrm{}..,u_r^{(n1)}(x),x)`$ (9)
$``$ $`\delta u_r^{(s)}(x)f_{u^{(s)}}u_(r1(x),u_{r1}^{(1)}(x),\mathrm{}..,u_{r1}^{(n1)}(x),x)].`$ (10)
The boundary conditions are similarly given by the difference of Eqs. 6 and 7 for two subsequent iterations:

$`{\displaystyle \underset{s=0}{\overset{n1}{}}}[\delta u_{r+1}^{(s)}(0)g_{ku^{(s)}}(u_r(0),u_r^{(1)}(0),\mathrm{}..,u_r^{(n1)}(0),0)`$ (11)
$`\delta u_r^{(s)}(0)g_{ku^{(s)}}(u_{r1}(0),u_{r1}^{(1)}(0),\mathrm{}..,u_{r1}^{(n1)}(0),0)]`$ $`=`$ $`0,`$ (12)
$`k`$ $`=`$ $`1,\mathrm{}l`$ (13)
and

$`{\displaystyle \underset{s=0}{\overset{n1}{}}}[\delta u_{r+1}^{(s)}(b)g_{ku^{(s)}}(u_r(b),u_r^{(1)}(b),\mathrm{}..,u_r^{(n1)}(b),b)`$ (14)
$`\delta u_r^{(s)}(b)g_{ku^{(s)}}(u_{r1}(b),u_{r1}^{(1)}(b),\mathrm{}..,u_{r1}^{(n1)}(b),b)]`$ $`=`$ $`0,`$ (15)
$`k`$ $`=`$ $`l+1,\mathrm{}n.`$ (16)
In view of the mean value theorem

$`f(u_r(x),u_r^{(1)}(x),\mathrm{}..,u_r^{(n1)}(x),x)f(u_{r1}(x),u_{r1}^{(1)}(x),\mathrm{}..,u_{r1}^{(n1)}(x),x)=`$ (17)
$`{\displaystyle \underset{s=0}{\overset{n1}{}}}\delta u_r^{(s)}(x)f_{u^{(s)}}(u_{r1}(x),u_{r1}^{(1)}(x),\mathrm{}..,u_{r1}^{(n1)}(x),x)+`$ (18)
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{s,t=0}{\overset{n1}{}}}\delta u_r^{(s)}(x)\delta u_r^{(t)}(x)f_{u^{(s)}u^{(t)}}(\overline{u}_{r1}(x),\overline{u}_{r1}^{(1)}(x),\mathrm{}..,\overline{u}_{r1}^{(n1)}(x),x),`$ (19)
where $`\overline{u}_{r1}^{(s)}(x)`$ lies between $`u_r^{(s)}(x)`$ and $`u_{r1}^{(s)}(x)`$. Now Eq. 10 can be written as

$`L^{(n)}\delta u_{r+1}(x){\displaystyle \underset{s=0}{\overset{n1}{}}}\delta u_{r+1}^{(s)}(x)f_{u^{(s)}}(u_r(x),u_r^{(1)}(x),\mathrm{}..,u_r^{(n1)}(x),x)=`$ (20)
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{s,t=0}{\overset{n1}{}}}\delta u_r^{(s)}(x)\delta u_r^{(t)}(x)f_{u^{(s)}u^{(t)}}(\overline{u}_{r1}(x),\overline{u}_{r1}^{(1)}(x),\mathrm{}..,\overline{u}_{r1}^{(n1)}(x),x).`$ (21)
Denoting $`G_r^{(n)}(x,y)`$ as the Greens function, which is the inverse of the following differential operator and incorporates linearized boundary conditions 6 and 7,

$`\stackrel{~}{L}^{(n)}=L^{(n)}{\displaystyle \underset{s=0}{\overset{n1}{}}}f_{u^{(s)}}(u_r(x),u_r^{(1)}(x),\mathrm{}..,u_r^{(n1)}(x),x){\displaystyle \frac{d^s}{dx^s}},`$ (22)
one can express the solution for the difference function $`\delta u_{r+1}`$ as

$`\delta u_{r+1}(x)=`$ (23)
$`{\displaystyle \frac{1}{2}}{\displaystyle _0^b}G_r^{(n)}(x,y)`$ $`{\displaystyle \underset{s,t=0}{\overset{n1}{}}}\delta u_r^{(s)}(y)\delta u_r^{(t)}(y)f_{u^{(s)}u^{(t)}}(\overline{u}_{r1}(y),\overline{u}_{r1}^{(1)}(y),\mathrm{}..,\overline{u}_{r1}^{(n1)}(y),y)dy.`$ (24)
The functions $`\delta u_r^{(s)}(y)\delta u_r^{(t)}(y)`$ could be taken outside of the sign of the integral at some point $`y=\overline{x}`$ belonging to the interval, so one obtains

$`\delta u_{r+1}(x)={\displaystyle \frac{1}{2}}{\displaystyle \underset{s,t=0}{\overset{n1}{}}}\delta u_r^{(s)}(\overline{x})\delta u_r^{(t)}(\overline{x})M_{st}(x).`$ (25)
where $`M_{sr}(x)`$ equals

$`M_{st}(x)={\displaystyle _0^b}G_r^{(n)}(x,y)f_{u^{(s)}u^{(t)}}(\overline{u}_{r1}(y),\overline{u}_{r1}^{(1)}(y),\mathrm{}..,\overline{u}_{r1}^{(n1)}(y),y)dy`$ (26)
If $`M_{st}(x)`$ is a strictly positive (negative) matrix for all $`x`$ in the interval, then $`\delta u_{r+1}(x)`$ will be positive (negative), and the monotonic convergence from below (above) results.

Obviously, from Eq. 24 follows

$`|\delta u_{r+1}(x)|k_r(x)\delta u_r^\mathit{2}`$ (27)
where $`k_r`$ is given by

$`k_r(x)={\displaystyle \frac{\mathit{1}}{\mathit{2}}}`$ $`{\displaystyle _0^b}|G_r^{(n)}(x,y)|{\displaystyle \underset{s,t=0}{\overset{n1}{}}}|f_{u^{(s)}u^{(t)}}(\overline{u}_{r1}(y),\overline{u}_{r1}^{(1)}(y),\mathrm{}..,\overline{u}_{r1}^{(n1)}(y),y)|dy`$ (28)
and $`\delta u_r`$ is a maximal value of any of $`|\delta \overline{u}_r^{(s)}|`$ on the interval (0,b).

Since Eq. 27 is correct for any $`x`$ on the interval (0,b), it is correct also for some $`x=\overline{x}`$ where $`|\delta u_{r+1}(x)|`$ reaches its maximum value $`\delta u_{r+1}`$. One therefore has

$`||\delta u_{r+1}|||k_r(\overline{x})||\delta u_r||^\mathit{2}`$ (29)
Assuming the boundness of the integrand in expression 28, that is the existence of the bounding function $`F(x)`$ such that integrand at $`x=\overline{x}`$ and at any $`y`$ is less or equal to $`F(y)`$, one finally has

$`||\delta u_{r+1}|||k||\delta u_r||^\mathit{2},`$ (30)
where

$`k={\displaystyle _\mathit{0}^b}F(x)\mathrm{𝑑𝑥}.`$ (31)
The linearized boundary conditions 6 and 7 are obtained from exact boundary conditions 2 and 3 by using the mean value theorem Eq. 19 and neglecting the quadratic terms, so that the error in using linearized boundary conditions vis-a-vis the exact ones is quadratic in the difference between the exact and linearized solutions. The maximum difference between boundary conditions 6 and 7 corresponding to two subsequent quasilinear iterations is therefore quadratic in $`\delta u_r`$. In view of this result and of Eq. 30, the difference between the subsequent iterative solutions of Eq.5 with boundary conditions 6 and 7 decreases quadratically with each iteration. In a similar way, one can show  that the difference $`\mathrm{\Delta }u_{r+1}(x)=u(x)u_r(x)`$ between the exact solution and the $`r`$-th iteration is decreasing quadratically as well:

$`||\mathrm{\Delta }u_{r+1}|||k||\Delta u_r||^\mathit{2}.`$ (32)
A simple induction of Eq. 30 shows  that $`\delta u_{n+1}(x)`$ for an arbitrary $`l<r`$ satisfies the inequality

$$\delta u_{r+1}(k\delta u_{l+\mathit{1}})^{\mathit{2}^{rl}}/k,$$
(33)
or, for $`l=0`$, we can relate the $`n+1`$ th order result to the 1st iterate by

$$\delta u_{n+1}((k\delta u_\mathit{1})^{\mathit{2}^n}/k.$$
(34)
The convergence depends therefore on the quantity $`q_1=ku_\mathit{1}u_\mathit{0}`$, where, as we have mentioned earlier, the zeroth iteration $`u_0(x)`$ is chosen from physical and mathematical considerations. Usually it is advantageous (see discussion below) that $`u_0(x)`$ would satisfy at least one of the boundary conditions. From Eq. (33) it follows, however, that for convergence it is sufficient that just one of the quantities $`q_m=k\delta u_m`$ is small enough. Consequently, one can always hope  that even if the first convergent coefficient $`q_1`$ is large, a well chosen initial approximation $`u_0`$ results in the smallness of at least one of the convergence coefficients $`q_m,m>1`$, which then enables a rapid convergence of the iteration series for $`r>m`$. It is important to stress that in view of the quadratic convergence of the QLM method, the difference $`\mathrm{\Delta }u_{r+1}`$ between the exact solution and the QLM iteration always converges to zero if the difference $`\delta u_{r+1}(x)`$ between two subsequent QLM iterations becomes infinitesimally small.

Indeed, if $`\delta u_r(x)`$ is close to zero, it means, since $`\delta u_{r+1}(x)=\mathrm{\Delta }u_r(x)\mathrm{\Delta }u_{r+1}(x)`$ that $`\mathrm{\Delta }u_r(x)=\mathrm{\Delta }u_{r+1}(x)`$ or $`Q_r=Q_{r+1}`$ where $`Q_r=k\Delta u_r`$. When one assumes the possibility that $`Q_r`$ and $`Q_{r+1}`$ could be not small, one could conclude that the iteration process “stagnates”, which means convergence to the wrong answer or no convergence at all.

However, such a conclusion is wrong since Eq. 32, which can be written as $`Q_{r+1}Q_r^2`$, for $`Q_r1`$ (this last inequality, starting from some r is a necessary condition of the convergence) could be not satisfied unless both $`Q_{r+1}`$ and $`Q_r`$ equal to zero. This proves that stagnation of the iteration process is impossible and convergence of $`\delta u_{r+1}`$ to zero automatically leads to convergence of the QLM iteration sequence to the exact solution. Hence the QLM assures not only convergence, but also convergence to the correct solution.

Another corollary of this iteration process is that if the solution and its derivatives are continuous functions of $`x`$, the convergence of the QLM in the whole region will follow. Indeed, even if the zero iteration $`u_0(x)`$ is chosen not to satisfy the boundary conditions, the next iteration $`u_1(x)`$, being a solution of a linear equation with linearized boundary conditions 6 and 7, will automatically satisfy the exact boundary conditions 2 and 3, at least up to the second order in difference $`\delta u_1`$ at the boundaries. This means that the difference between the exact and first QLM iterations at some intervals near the boundaries will be small, so that the QLM iterations in this interval would converge. Because the subsequent values of $`k\delta u_m(x),m>\mathit{2}`$ became much smaller for this interval, in view of assumed continuity of the solution and its derivatives these differences will also be small at the neighboring intervals. The subsequent iterations will extend the convergence to the next neighboring intervals and so on, until the convergence in the whole region will be reached. The predicted trend is therefore that the QLM yields rapid convergence starting at the regions where the boundary conditions are imposed and then spreading from there to all other regions.

An additional important corollary is that, in view of Eq. 33, once the quasilinear iteration sequence starts to converge, it will continue to do so, unlike the perturbation expansion, which is often given by an asymptotic series and therefore converges only up to a certain order and diverges thereafter.

Based on this summary of the QLM, one can deduce the following important features of the quasilinearization method:

* The method approximates the solution of nonlinear differential equations by treating the nonlinear terms as a perturbation about the linear ones, and is not based, unlike perturbation theories, on the existence of some kind of small parameter.
* The iterations converge uniformly and quadratically to the exact solution. In case of matrix $`M_{st}`$ in Eq. 26 being strictly positive (negative) for all $`x`$ in the interval, the convergence is also monotonic from below (above).
* For rapid convergence it suffices that an initial guess for the zeroth iteration is sufficiently good to ensure the smallness of just one of the quantities $`q_r=ku_{r+\mathit{1}}u_r`$. If the solution and its derivatives are continuous, convergence follows from the fact that starting from the first iteration, all QLM iterations automatically satisfy the quasilinearized boundary conditions 6 and 7. The convergence is extremely fast: if, for example, $`q_1`$ is of the order of $`\frac{1}{3}`$, only 4 iterations are necessary to reach the accuracy of 8 digits, since $`(\frac{1}{3})^{2^n}`$ is of the order of $`(\frac{1}{10})^{2^{n1}}`$.
* Convergence of $`\delta u_{r+1}`$ to zero automatically leads to convergence of the QLM iteration sequence to the exact solution.
* Once the quasilinear iteration sequence at some interval starts to converge, it will always continue to do so. Unlike an asymptotic perturbation series, the quasilinearization method yield the required precision once a successful initial guess generates convergence after a few steps.

## III Analytically solvable example: comparison of quasilinearization approach with exact solution and with perturbation theory.

In order to investigate the applicability of the quasilinearization method and its convergence and numerical stability, let us start from a simple example of an analytically solvable nonlinear ordinary differential equation suggested in ref. :

$$u^{}(r)=gu^n(r),u(0)=1,$$
(35)
where the boundary condition at $`r=0`$ is also given and  means differentiation in variable $`r`$. The exact solution to this problem is

$$u(r)=(1+(n1)gr)^{\frac{1}{n1}}$$
(36)
Since

$$(1+x)^q=\underset{0}{\overset{\mathrm{}}{}}\frac{\mathrm{\Gamma }(q+1)}{m!\mathrm{\Gamma }(q+1m)}x^m,$$
(37)
the expansion of the solution 36 in powers of $`g`$ is given by

$$u(r)=\underset{0}{\overset{\mathrm{}}{}}\frac{\mathrm{\Gamma }(\frac{n2}{n1})(g(n1))^m}{m!\mathrm{\Gamma }(\frac{n2}{n1}m)}r^m$$
(38)
The convergence radius of the series 38 is $`R=1/(g(n1)),`$ which is inversely proportional to the extent $`n1`$ of the nonlinearity and to the value $`g`$ of the perturbation parameter.

Now consider the quasilinearization approach to this equation, taking, for example, $`g=1`$ and $`n=6.`$ Here we consider Eq. 35 with a rather strong degree of nonlinearity. In this case, one can expect the convergence of the perturbation expansion only up to $`rR=\frac{1}{5}`$.

The QLM procedure in the case where the nonlinear term depends only on the solution itself and not on its derivatives reduces to setting $`u_{k+1}^{}(r)=f(u_k)+(u_{k+1}(r)u_k(r))f_u(u_k)`$. Here $`f=gu^n(r)`$ while its functional derivative $`f_u`$ equals to $`gnu^{n1}(r)`$. The quasilinearized equation 5 for the $`(k+1)`$-th iteration for this case has therefore the following form:

$$u_{k+1}^{}(r)+ngu_k^{n1}(r)u_{k+1}(r)=(n1)gu_k^n(r),u_{k+1}(0)=1,$$
(39)
where $`u_k(r)`$ is a previous iteration which is considered to be a known function. Let us choose as a zero iteration $`u_0(r)1`$ which satisfies the boundary condition $`u_0(0)=1`$.

The results of our QLM calculations with Eq. 39 are presented in Fig. 1 which displays the exact solution for the case of $`n=6`$ and $`g=1,`$ together with the first four QLM iterations. Convergence to the exact solution in Fig. 1 is monotonic from above as it should be as discussed in Section II and in Refs. due to fact that the second functional derivative $`n(n1)u^{n2}(x)`$ of the left-hand side of Eq.35 for even $`n`$ is strictly negative. The convergence starts at the boundary, exactly as expected from the discussion in section II, and expands with each iteration to larger values of the variable $`r`$. The difference between the exact solution and the sixth QLM iteration for all r in the range between zero and five where our calculations were performed is less than $`10^6.`$ Note that the QLM yields a solution beyond the convergence radius limit on the series solution of $`1/5.`$

## IV Lane-Emden equation

The Lane-Emden equation

$$y^{\prime \prime }(r)+\frac{2}{r}y^{}(r)+y^n(r)=0,y(0)=1,y^{}(0)=0$$
(40)
describes a variety of phenomena in theoretical physics and astrophysics, including aspects of stellar structure, the thermal history of a spherical cloud of gas, isothermal gas spheres, and thermionic currents, see ref.  and references therein. The parameter $`n`$ defines an equation of state with its physically interesting range being $`0n5`$. The equation also appears in other contexts, e.g., in case of radiatively cooling, self gravitating gas clouds, in the mean-field treatment of a phase transition in critical absorption or in the modeling of clusters of galaxies. For $`n=0,1`$ and $`5`$ the equation can be solved analytically. Setting $`y=\frac{u}{r}`$ transforms the equation to a more convenient form without a first derivative:

$$u^{\prime \prime }(r)+\frac{u^n(r)}{r^{n1}}=0,u(0)=0,u^{}(0)=1.$$
(41)
Let us consider this nonlinear equation for the physically interesting and analytically nonsolvable case of $`n=4`$. The quasilinearized form of equation 41 is

$$u_{k+1}^{\prime \prime }(r)+n\frac{u_k^{n1}(r)}{r^{n1}}u_{k+1}(r)=\frac{n1}{r^{n1}}u_k^n(r),u_{k+1}(0)=0,u_{k+1}^{}(0)=1.$$
(42)
The simplest initial guess, satisfying the boundary conditions will be $`u_0(r)=r.`$ Comparison of the quasilinear solutions corresponding to the first five iterations with the numerically computed exact solution are given in Fig. 2. The figure shows that the convergence to the exact solution is very fast. It starts, as in the example of the previous section, at the left boundary and spreads with each iteration to larger values of $`r`$ as expected from the discussion in section II. The difference between the exact solution and the eighth QLM iteration for all $`r`$ in the range between zero and ten, where our calculations were performed, is less than $`10^{11}`$ .

## V Thomas-Fermi equation

The Thomas-Fermi equation

$$\sqrt{x}u^{\prime \prime }(x)=u^{\frac{3}{2}}(x),u(0)=1,u(\mathrm{})=0,$$
(43)
is an equation for the electron density around the nucleus of the atom. The left hand side of the above equation equals zero for $`u<0.`$ The Thomas-Fermi equation is also very useful for calculating form-factors and for obtaining effective potentials which can be used as initial trial potentials in self-consistent field calculations.. It is also applicable to the study of nucleons in nuclei and electrons in metal. It is long known (see  and references therein) that solution of this equation is very sensitive to a value of the first derivative at zero which insures smooth and monotonic decay from $`u(0)=1`$ to $`u(\mathrm{})=0`$ as demanded by boundary conditions. By contrast, the computation is much simpler for the quasilinearized version of this equation. The QLM procedure in this case reduces to setting $`u_{k+1}^{\prime \prime }(r)=f(u_k)+(u_{k+1}(r)u_k(r))f_u(u_k)`$, where $`f=\frac{u^{3/2}(r)}{\sqrt{x}}`$ and the functional derivative is $`f_u=(3/2)\frac{u^{1/2}}{\sqrt{x}}`$, so that the QLM equation has a form:

$$\sqrt{x}u_{k+1}^{\prime \prime }(x)\frac{3}{2}u_k^{\frac{1}{2}}(x)u_{k+1}(x)=\frac{1}{2}u_k^{\frac{3}{2}}(x),u_{k+1}(0)=1,u_{k+1}(\mathrm{})=0,$$
(44)
which is easily solved by specifying directly the boundary condition at infinity without searching first for the proper value of the first derivative. The initial guess, satisfying the boundary condition at zero was chosen to be $`u_0(x)1`$. The results of QLM calculations with Eq. 44 are presented in Fig. 3 which displays the exact solution together with the first four QLM iterations. The convergence starts at the boundaries, exactly as expected from the discussion in section II, and expands with each iteration to a wider range of values of the variable $`x`$. The difference between the exact solution and the eighth QLM iteration for all $`x`$ in the range between zero and forty where our calculations were performed is less than $`10^7.`$

## VI Classical anharmonic oscillator

The classical anharmonic oscillator satisfies the nonlinear second-order Duffing equation

$$\ddot{u}(t)+u(t)+gu^3(t)=0.$$
(45)
The typical initial conditions are

$$u(0)=1,\dot{u}(0)=0.$$
(46)
The solution oscillates strongly and thus is more difficult to approximate. It is, for example, well known  that the usual perturbative solution is valid only for times $`t`$ small compared with $`\frac{1}{g}`$, so that for larger $`g`$ the perturbative solution is adequate only on a small time interval. In contrast, the quasilinearization approach gives solution in the whole region also for large $`g`$-values.

The quasilinearized equation is

$$\ddot{u}_{k+1}(t)+(1+3gu_k^2(t))u_{k+1}(t)2gu_k^3(t)=0,u_{k+1}(0)=1,\dot{u}_{k+1}(0)=0.$$
(47)
The results of QLM calculations with Eq. 47 for $`g=3`$ are presented in Figs. 4 and 5. Fig. 4 displays the exact solution together with the QLM solutions for the first, second and fourth iterations while Fig. 5 shows comparison of exact solution with sixth, seventh and eighth QLM iterations. Again, the convergence starts at the left boundary as expected from the discussion in section II, and expands with each iteration to larger values of the variable $`t`$. The difference between the exact solution and the eleventh QLM iteration for all t in the range between zero and seven where our calculations were performed is less than $`10^{10}`$ .

## VII Blasius equation

A third order nonlinear Blasius equation

$$u^{\prime \prime \prime }(x)+u^{\prime \prime }(x)u(x)=0,u(0)=u^{}(0)=0,u^{}(\mathrm{})=1$$
(48)
describes the velocity profile of the fluid in the boundary layer. The QLM procedure in this case is given by $`u_{k+1}^{\prime \prime \prime }(x)=f(u_k,u_k^{\prime \prime })+(u_{k+1}u_k)f_u(u_k,u_k^{\prime \prime })+(u_{k+1}^{\prime \prime }u_k^{\prime \prime })f_{u^{\prime \prime }}(u_k,u_k^{\prime \prime })`$, where $`f(u,u^{})=u^{\prime \prime }u`$, $`f_u(u,u^{})=u^{\prime \prime }`$ and $`f_{u^{\prime \prime }}(u,u^{})=u`$. The quasilinearized version of the Blasius equation thus has a form

$`u_{k+1}^{\prime \prime \prime }(x)+u_k(x)u_{k+1}^{\prime \prime }(x)+u_{k+1}(x)u_k^{\prime \prime }(x)u_k(x)u_k^{\prime \prime }(x)=0,`$ (49)
$`u_{k+1}(0)=u_{k+1}^{}(0)=0,u_{k+1}^{}(\mathrm{})=1.`$ (50)
The initial guess, satisfying the boundary condition for the derivative at zero was chosen to be $`u_0(x)1`$. The results of QLM calculations with Eq. 50 are presented in Fig. 6 which displays the exact solution together with the first QLM iteration. The convergence starts at the left boundary as follows from the discussion in section II, and expands with each iteration to larger values of the variable $`x`$. The difference between the exact solution and the fifth QLM iteration for all $`x`$ in the range between zero and ten where our calculations were performed is less than $`10^{11}`$ .

## VIII Conclusion

Summing up, we formulated here the conditions under which the quadratic, uniform and often monotonic convergence of the quasilinearization method are valid.

We have followed here the derivation outlined in ref. , which is not based, unlike the derivations in refs. , on a smallness of the interval and on the boundness of the nonlinear term and its functional derivatives, the conditions which usually are not fulfilled in physical applications.

In order to analyze and highlight the power and features of the quasilinearization method (QLM), in this work we have also made numerical computations on different ordinary second and third order nonlinear differential equations, applied in physics, such as the Blasius, Duffing, Lane-Emden and Thomas-Fermi equations and have compared the results obtained by the quasilinearization method with the exact solutions. Although all our examples deal only with linear boundary conditions, the nonlinear boundary conditions can be handled readily after their quasilinearization as explained in Section II.

In all the examples considered in the paper the simplest initial guess was enough to produce a quadratic convergence. However, as it is well known for the Newton method on which quasilinearization method is based that the divergence could occur if the initial guess is too bad. In this case the modification of QLM based on the damped Newton method with relaxation factor /cite CB,RR may help. The price of such modification would be, like in the damped Newton method only linear convergence instead of the quadratic one.

Our conclusions are as follows:

The QLM treats nonlinear terms by a series of nonperturbative iterations and is not based on the existence of some kind of small parameter. At every iterative stage, the differential operator changes significantly to account for the nonlinearity, which is the major way that the QLM differs from other approximative techniques. As a result, as we see in all our examples, the QLM is able to handle large values of the coupling constant and any degree of the nonlinearity, unlike perturbation theory. Thus the QLM provides extremely accurate and numerically stable answers for a wide range of nonlinear physics problems. The QLM is also very easy to apply.

In view of all this, since most equations of physics, from classical mechanics to quantum field theory, are either not linear or could be transformed into a nonlinear form, the quasilinearization method appears to be extremely useful and in many cases more advantageous than the perturbation theory or its different modifications, like expansion in inverse powers of the coupling constant, the $`1/N`$ expansion, etc.

###### Acknowledgements.

The authors are indebted to Dr. Schönauer for careful reading of the manuscript and many constructive comments and suggestions. The research was supported in part by the U.S. National Science Foundation grant PHY-9970775 (FT) and by the Israel Science Foundation grant 131/00 (VBM).

FIGURE CAPTIONS

FIG. 1. Convergence of QLM iterations for the analytic example of section III and comparison with the exact solution. Thin solid, dot-dashed, short-dashed and dotted curves correspond to the first, second, third and fourth QLM iteration respectively, while the thick solid curve displays the exact solution. The convergence is monotonic from above as it should be according to the discussion in the text. The difference between the exact solution and the sixth QLM iteration for all r in the figure is less than $`10^6`$.

FIG. 2. Convergence of QLM iterations for the Lane-Emden equation and comparison with the numerically obtained exact solution. Thin solid, dot-dashed, short-dashed, long-dashed and dotted curves correspond to the first, second, third, fourth and fifth QLM iteration, respectively, while the thick solid curve displays the exact solution. The difference between the exact solution and the eighth QLM iteration for all r in the figure is less than $`10^{11}`$

FIG. 3. Convergence of QLM iterations for the Thomas-Fermi equation and comparison with the numerically obtained exact solution. Thin solid, dot-dashed, short-dashed and dotted curves correspond to the first, second, third and fourth QLM iteration, respectively, while the thick solid curve displays the exact solution. The difference between the exact solution and the eighth QLM iteration for all x in the figure is less than $`10^7`$ .

FIG. 4. Convergence of the first few QLM iterations for the Duffing equation and comparison with the numerically obtained exact solution. The dotted curves on three consecutive graphs correspond to the first, second and fourth QLM iteration respectively, while the solid curve displays the exact solution.

FIG. 5. Convergence of the higher QLM iterations for the Duffing equation and comparison with the numerically obtained exact solution. The dotted curves on the three consecutive graphs correspond to the sixth, seventh and eighth QLM iteration respectively, while the solid curve displays the exact solution. The difference between the exact solution and the eighth QLM iteration for all t in the figure is less than $`10^{10}`$ .

FIG. 6. Comparison of the first QLM iteration for the Blasius equation with the numerically obtained exact solution. The difference between the exact solution and the fifth QLM iteration for all x in the figure is less than $`10^{10}`$ .


text:
# Heavy-Fermion Formation at the Metal-to-Insulator Transition in Gd1-xSrxTiO3

## Abstract

The perovskite-like transition-metal oxide Gd1-xSrxTiO3 is investigated by measurements of resistivity, specific-heat, and electron paramagnetic resonance (EPR). Approaching the metal-to-insulator transition from the metallic regime ($`x0.2`$), the Sommerfeld coefficient $`\gamma `$ of the specific heat becomes strongly enhanced and the resistivity increases quadratically at low temperatures, which both are fingerprints of strong electronic correlations. The temperature dependence of the dynamic susceptibility, as determined from the Gd3+-EPR linewidth, signals the importance of strong spin fluctuations, as observed in heavy-fermion compounds.

PACS: 71.30, 71.27, 76.30

The transition-metal oxides $`R`$TiO3 (where $`R`$ denotes Y3+, La3+, or some trivalent rare-earth ion) are known as typical Mott-Hubbard insulators with Ti3+ in a 3$`d^1`$ electronic configuration . Doping $`R_{1x}A_x`$TiO3 with divalent alkali-earth ions $`A`$ like Ca2+ or Sr2+ partially changes the Ti valence to Ti4+ (3$`d^0`$) and induces metallic behavior at a certain doping concentration $`x`$. The metal-to-insulator transition of La1-xSrxTiO3 and Y1-xCaxTiO3 has been investigated in detail . A critical increase of the Sommerfeld coefficient $`\gamma `$ of the specific heat approaching the metal-to-insulator transition has been observed in both compounds with $`\gamma 25`$ mJ/(mole K2) in Y0.6Ca0.4TiO3, indicating an enhancement of the effective electron mass. Heavy-fermion formation has been observed in a number of transition-metal oxides, with Sr1-xCaxRuO3  and LiV2O4  being the most prominent examples. In the latter two compounds the temperature dependence of the dynamic susceptibility, as probed by the spin-lattice relaxation time in NMR experiments, played a key role in identifying the heavy-fermion behavior. In this letter we investigate the dynamic susceptibility utilizing Gd EPR and provide experimental evidence of heavy-fermion behavior in Gd1-xSrxTiO3.

Ceramic samples of Gd1-xSrxTiO3 have been calcinated from TiO2, SrCO3, Gd2O3 and Ti2O3 powders of high purity (better than 3N) and pressed into pellets. These pellets have been annealed at 1473 K under N2 atmosphere for 15 hours and finally arc-melted under argon atmosphere. For Sr concentrations $`x0.7`$ X-ray powder diffraction revealed the proper orthorhombically distorted GdFeO3 perovskite structure . For $`x>0.7`$ the system attains the cubic perovskite structure like pure SrTiO3. The volume of the unit cell remains nearly constant at about 235 $`\AA ^3`$ in the whole concentration range. Susceptibility measurements were performed with a commercial SQUID magnetometer (Quantum Design) in a temperature range 4.2 $`T`$ 400 K. Pure GdTiO3 shows a ferrimagnetic susceptibility, which follows a Curie-Weiss law $`\chi (T+\mathrm{\Theta })^1`$ with $`\mathrm{\Theta }=10`$ K at high temperatures $`T>100`$ K, and exhibits an ordering temperature $`T_C28`$ K. This fits reasonably to the value of 34 K given in literature , where the ordered state is reported to show an antiparallel alignment of the ferromagnetically ordered Ti3+ spins with respect to the Gd3+ spins. With increasing $`x`$ the Curie-Weiss temperature $`\mathrm{\Theta }`$ shows approximately a linear decrease and vanishes for $`x1`$, whereas the magnetic order strongly becomes suppressed. Already at $`x>0.1`$ no magnetic order was observed for $`T1.7`$ K from our EPR measurements reported below. For all samples, the paramagnetic susceptibility is dominated by the contribution of the Gd3+ spins ($`S_{\mathrm{Gd}}`$ = 7/2), which is by a factor $`S_{\mathrm{Gd}}(S_{\mathrm{Gd}}+1)/S_{\mathrm{Ti}}(S_{\mathrm{Ti}}+1)`$ = 21 times stronger than the contribution of the Ti3+ spins ($`S_{\mathrm{Ti}}`$ = 1/2). Comparing the theoretically expected and experimentally observed susceptibility, we obtained that the samples exhibit the correct composition $`x`$ within an error of $`\mathrm{\Delta }x0.02`$.

Electrical resistivity, using standard four-probe lock-in technique, and specific heat, using the adiabatic method, were measured in home-built 4He cryostats in a temperature range 1.5 $`T`$ 300 K on bulk samples cut from the arc-melted ingots. Fig. 1 shows the temperature dependence of the resistivity for Sr concentrations $`x0.5`$ normalized to the value at $`T=300`$ K. For pure GdTiO3 the resistivity strongly increases on decreasing temperature, as it is expected for a Mott insulator. With increasing Sr concentration the resistivity strongly decreases and a metal-to-insulator transition occurs for $`x0.2`$. Within the concentration range 0.3 $`x`$ 0.5 the resistivity obeys roughly a quadratic temperature dependence, which is demonstrated in the inset of Fig. 1. The observed temperature dependence $`\rho \rho _0T^2`$, instead of the usual phonon contribution $`\rho \rho _0T^5`$, is typical for electron-electron scattering and only observable, if the effective mass of the electrons is strongly enhanced due to electronic correlations like in heavy-fermion systems . At Sr concentrations $`x>0.6`$ (not shown in Fig. 1) the system becomes gradually insulating, approaching the band insulator SrTiO3 .

Fig. 2 shows the temperature dependence of the specific heat $`C`$ for selected samples. Using the representation $`C/T`$ as a function of $`T^2`$, the linear increase is due to the phonon contribution $`\beta `$, following $`C_{\mathrm{ph}}=\beta T^3`$, and its extrapolation to $`T=0`$ yields the Sommerfeld coefficient $`\gamma `$ due to the electronic contribution $`C_{\mathrm{el}}=\gamma T`$. The increase of the data below 10 K is due to the onset of magnetic order of the Gd spins below 1 K, as we proved for $`x=0.4`$ in a 3He/4He-dilution refrigerator. From the extrapolation of the data above 10 K down to $`T=0`$ we determined the Sommerfeld coefficient $`\gamma `$ for the metallic regime, which is shown in the inset of Fig. 2. Approaching the metal-to-insulator transition from the metallic side, the Sommerfeld coefficient is strongly enhanced up to 50 times with respect to simple metals at $`x=0.2`$. This is another hint concerning the enhancement of the effective electronic mass $`m^{}`$ or - in the same sense - of the density of states $`N(E_\mathrm{F})`$ at the fermi energy, because $`\gamma m^{}N(E_F)`$. Comparing these results to resistivity and specific heat in the well investigated systems La1-xSrxTiO3 and Y1-xCaxTiO3, we find a quite good agreement. An enhancement of the electronic density of states at low temperatures around the Mott-Hubbard transition is also predicted by theoretical calculations . Hence Gd1-xSrxTiO3 is an appropriate system to study the dynamic susceptibility close to a

Mott-Hubbard transition.

The EPR experiments have been performed at a Bruker ELEXSYS E500-CW spectrometer at X-band frequencies (9.4 GHz) equipped with continuous gas-flow cryostats for He (Oxford Instruments) and N2 (Bruker) in the temperature range between 4.2 and 680 K. For temperatures down to 1.7 K we used a cold-finger 4He-bath cryostat. To avoid the influence of the skin effect, which prevents the microwave to enter metallic samples, the arc-melted ingots were powdered, filled in quartz tubes, and embedded in paraffin or NaCl for temperatures below or above room temperature, respectively.

The EPR spectra of all compounds consist of a single exchange-narrowed resonance line of Lorentzian shape with a g value $`g2`$. The evaluation of comparable Gd3+ spectra has been described in detail in a previous publication . Their intensity follows the Curie-Weiss behavior of the static susceptibility measured with the SQUID magnetometer and comparison with a standard sample showed that it is due to all Gd3+ spins. In the Mott-insulating regime $`x0.15`$ the resonance linewidth $`\mathrm{\Delta }H`$ strongly increases with decreasing temperature below $`T<50`$ K, due to the onset of magnetic order. In the minimum near 50 K, it attains a value of about 1 kOe. To higher temperatures $`T>50`$ K the linewidth increases with a slightly negative curvature and approximately 2.5 Oe/K. In the metallic regime $`0.2x0.6`$, the linewidth $`\mathrm{\Delta }H`$ exhibits a clearly different behavior, which is shown in Fig. 3: Starting at low temperatures the linewidth strongly increases, develops a prominent maximum and decreases again. At temperatures above 150 K the linewidth remains nearly constant ($`x0.4`$) or shows a further slight decrease ($`x0.5`$). Approaching the metal-to-insulator transition from above ($`x0.2`$), the maximum increases and shifts to lower temperatures.

In the following discussion, we confine ourselves to the metallic regime: In metals the transversal spin-relaxation time $`T_2`$ equals the longitudinal or spin-lattice relaxation time $`T_1`$. Hence the EPR linewidth $`\mathrm{\Delta }H1/T_2`$ directly measures the spin-lattice relaxation rate $`1/T_1`$, which is usually found linear in temperature for simple metals, probing the electronic density of states $`N(E_\mathrm{F})`$ due to the Korringa law :

$$\mathrm{\Delta }H=bT,\mathrm{where}bN^2(E_\mathrm{F}).$$
(1)
However, the temperature evolution of the linewidth observed in Gd1-xSrxTiO3 differs completely from the linear Korringa law , but it reveals similar features as the Gd ESR  and NMR spin-lattice relaxation rate  in heavy-fermion compounds. Moreover comparable NMR results have been reported form transition-metal oxides with heavy-fermion formation like LiV2O4  and Ca1-xSrxRuO3 . In all these systems a second important contribution to the spin-lattice relaxation arises from spin fluctuations of the electronic bath due to the fluctuation-dissipation theorem, which for temperatures large compared to the exciting frequency $`k_BT\mathrm{}\omega `$ reads

$$\frac{1}{T_1}\frac{T}{\omega }\mathrm{Im}\chi (\omega ,T).$$
(2)
Here Im$`\chi (\omega ,T)`$ denotes the imaginary part of the dynamic susceptibility of the band states represented by the Ti $`3d`$ electrons in the present case. An appropriate expression for the dynamic susceptibility, which takes the itinerant character of the Ti electrons into account, is given by Ishigaki and Moriya , who developed a theoretical description of nuclear magnetic relaxation around the magnetic instabilities in metals. Their results were also used to explain the NMR relaxation rate in Ca1-xSrxRuO3 . Ishigaki’s simulations of the nuclear spin-lattice relaxation rate $`1/T_1`$ near a ferromagnetic instability - the Ti spins order ferromagnetically (!) - look quite similar to our linewidth data. Substituting the hyperfine coupling of the nuclear spin with the electronic system by the super-exchange interaction between Gd3+ and Ti3+ spins, we achieved an analogous expression for the EPR linewidth. Moreover it was necessary, to include the effect of the external magnetic field onto the magnetic instability . Finally we used the following expression to fit our data:

$$\mathrm{\Delta }H=A_{\mathrm{SE}}\frac{3ty^2}{2(y^3+h)}+\mathrm{\Delta }H_0.$$
(3)
The reduced inverse 3$`d`$ susceptibility $`y`$ is given by an implicit integral equation , which has to be iterated numerically. It is characterized by two parameters $`y_0`$ and $`y_1`$, respectively: $`y_0`$ is proportional to the inverse static susceptibility at zero temperature and therefore vanishes at the magnetic instability. $`y_1T_0y_2`$ is proportional to the the energy width $`T_0`$ ($`0.010.1`$ eV) of the dynamical spin-fluctuation spectrum. Both parameters are normalized to the width $`T_\mathrm{A}`$ ($`1`$ eV) of the distribution of the static susceptibility in momentum space. The variable $`t=T/T_0`$ denotes the reduced temperature. The parameter $`hH^2`$ describes the influence of the applied magnetic field $`H`$. The super-exchange coupling, which is included in the prefactor $`A_{\mathrm{SE}}`$, is assumed to be temperature independent. The residual linewidth $`\mathrm{\Delta }H_0`$ is due to the contribution of impurities and inhomogeneties in the sample. Hence there are six parameters altogether to describe the temperature dependence of the linewidth.

To obtain a reasonable description of our data, we tried to find a parameter set, where the parameters $`h`$ \- the resonance is found at about $`H=3.4`$ kOe for all compounds - and $`y_2(1/T_\mathrm{A})^2`$, which is related to the width of the conduction band, attain the same values for all Sr concentrations. The coupling constant $`A_{\mathrm{SE}}`$ and the residual linewidth $`\mathrm{\Delta }H_0`$ were allowed to change slightly. The most important effect was expected to appear in the parameter $`y_0`$, which is directly tuned by the Sr concentration $`x`$, and the energy width $`T_0`$ of the spin-fluctuation spectrum. The results are drawn as solid lines in Fig. 3 and the respective fit-parameters are listed in Table 1. The fit nicely describes the experimental data at temperatures $`TT_0`$ and underlines the dominant influence of the 3$`d`$-spin fluctuations for relaxation near the Mott-Hubbard transition: Indeed the super-exchange coupling and residual linewidth change only slightly within the metallic concentration regime under consideration. Both parameters $`y_0`$ and $`T_0`$ increase nearly linearly with the Sr concentration $`x`$. The critical concentration for the onset of magnetic order ($`y_0=0`$) accompanied by the transition to the Mott insulator is found at $`x0.2`$ in agreement with the resisitivity measurements.

The resonance linewidth at temperatures far below the maximum, which is shown in the upper inset of Fig. 3, exhibits deviations from the linear behavior predicted by the theory. The slope is too small, and for $`x0.4`$ one observes even a slight upturn with decreasing temperature instead of the expected linear decrease. These deviations are probably due to the influence of the crystal-electric field or the onset of magnetic order below $`T=1`$ K, which both induce a broadening of the resonance line. In the case of a small crystal-electric field the linewidth at low

temperatures follows (solid lines in the upper inset of Fig. 3)

$$\mathrm{\Delta }H=b_{\mathrm{LT}}T+\frac{M_2}{b_{\mathrm{LT}}T}+\mathrm{\Delta }H_0,$$
(4)
where $`M_2`$ is proportional to the second moment of the crystal-field splitting. It is interesting to compare the low-temperature Korringa slope $`b_{\mathrm{LT}}`$ derived from equation (4) - having in mind the pure Korringa law (1) - with the squared Sommerfeld coefficient $`\gamma ^2`$ of the specific heat, as is shown in the lower inset of Fig. 3. The clear correlation of both parameters underlines the enhancement of the electronic density of states $`N(E_\mathrm{F})`$ at low temperatures in the metallic regime near the Mott-Hubbard transition. It is visible by both macroscopic ($`\gamma `$) and microscopic ($`b_{\mathrm{LT}}`$) probes.

In conclusion we have shown that the compound Gd1-xSrxTiO3 exhibits a metal-to-insulator transition near a Sr concentration $`x=0.2`$: Electrical resistivity and specific heat at the transition are in full agreement with the well known systems La1-xSrxTiO3 and Y1-xCaxTiO3. As Gd3+ naturally belongs to the compound, it is well suited for EPR investigations. The temperature dependence of the Gd3+ linewidth in the metallic regime $`0.2x0.6`$ is determined by the Ti3+-spin fluctuations and is well described by the theory of Ishigaki and Moriya for a metal near a ferromagnetic instability. The low-temperature behavior of resistivity, specific heat and ESR linewidth resembles to the properties of heavy-fermion systems, but the enhancement of the effective electronic masses by a factor 50 is weaker than in f-derived heavy fermions, which are dominated by Kondo-compensation effects and where one observes an enhancement factor of 1000.

We are grateful to Dr. E. W. Scheidt, Experimentalphysik III, University of Augsburg, for specific-heat measurements at temperatures below 1 K. This work was supported by the Bundesministerium für Bildung und Forschung (BMBF) under the contract number EKM 13N6917/0 and by the Deutsche Forschungsgemeinschaft within Sonderforschungsbereich 484.