text:
# A Monte Carlo Test of the Fisher-Nakanishi-Scaling Theory for the Capillary Condensation Critical Point

## I Introduction

The application prospects of nanoscale technology have created a fresh interest in the behavior of both simple fluids and complex fluids confined in pores or in a thin film geometry in layers confined by parallel walls . However, a prerequisite for the clarification of pattern formation  and dynamics  is a good understanding of the interplay between bulk and surface effects on thermodynamics and phase behavior in this finite-size geometry : although theoretical aspects of phase transitions and critical phenomena in confined geometry have been considered since a long time , this is still a topic of active current research  even for one of the most well-known phenomena, namely “capillary condensation” . By “capillary condensation” one means the finding, already discovered in the 19th century , that in a capillary the condensation of a gas occurs already at a lower pressure $`p`$ than the coexistence pressure $`p_{\text{coex}}`$ necessary to induce condensation in the bulk. Qualitatively, this shift of the transition can be attributed to the interaction of the fluid molecules with the attractive walls of the capillary. Although confinement effects on fluids and their phase transitions have been studied experimentally since a long time as well , a quantitative characterization of the shift of the capillary condensation critical point remains a challenge. While for temperatures $`T`$ below the critical temperature $`T_c(D)`$ of the thin film of thickness $`D`$ the chemical potential at the condensation transition $`\mu _c(D)`$ is shifted relative to its bulk value simply as $`\mu _c(D)\mu _c(\mathrm{})D^1`$ (“Kelvin equation”) , for large enough $`D`$, Fisher and Nakanishi  predicted a completely different power law for the corresponding shift at the critical temperature itself, namely

$$\mu _c(D)\mu _c(\mathrm{})D^{(\mathrm{\Delta }\mathrm{\Delta }_1)/\nu },T=T_c(D)$$
(1)
for weak surface forces. In Eq. (1), critical exponents of the three-dimensional Ising model universality class (that encompasses criticality of gas-fluid critical points or the related unmixing transitions in binary mixtures, etc.) enter, namely the correlation length exponent  $`\nu 0.63`$ and the “gap exponent” $`\mathrm{\Delta }=\gamma +\beta 1.56`$ and the corresponding exponent for a free surface  $`\mathrm{\Delta }_10.460.48`$. Also for the shift of $`T_c`$ a similar power law holds,

$$T_c(\mathrm{})T_c(D)D^{1/\nu },$$
(2)
which is the same relation as is familiar from standard finite-size scaling  for the shift of $`T_c`$ in films with “neutral walls” {i.e., no surface field preferring one of the phases coexisting for $`T<T_c(D)`$ at $`\mu _c(D)\mu _c(\mathrm{})`$ act} or in films with periodic boundary conditions where surface effects are a priori absent. While for the latter systems Eq. (2) has been studied by various methods , Eq. (1) to our knowledge has not yet been tested by Monte Carlo simulations. In previous work , tests of the Kelvin equation and corrections to it  have been carried out and a capillary condensation critical point was located for a thin Ising film  but for a single value of $`D`$ only.

In the present paper, we fill in this gap by presenting a Monte Carlo study of the critical behavior of capillary condensation in thin Ising films for a range of thicknesses. Invoking the universality principle , one can argue that nearest neighbor Ising lattices with short range surface fields should yield the same power law, Eq. (1), as more realistic models and real fluids in slit-like capillaries do. Unlike the situation in real fluids, packing effects at the surfaces and a dependence of the density in the middle of the film on its thickness are, however, absent and one focuses on the universal critical behavior. One can argue that the order parameter correlations in the directions parallel to the wall should all scale with the critical exponents of the universality class of the two-dimensional Ising model ,

$`M\stackrel{~}{t}^{\beta _2},\chi |\stackrel{~}{t}|^{\gamma _2},\xi _{||}|\stackrel{~}{t}|^{\nu _2},\beta _2=1/8,\gamma _2={\displaystyle \frac{7}{4}},\nu _2=1,`$

$$\stackrel{~}{t}1T/T_c(D)0\text{,\hspace{1em}all }D<\mathrm{}.$$
(3)
Of course, one expects that for large $`D`$ the asymptotic critical region where this two-dimensional critical behavior holds is very narrow, due to a crossover to the three-dimensional critical behavior as $`D\mathrm{}`$, and a quantitative understanding of this crossover  is a challenging aspect of this problem, too.

In Sec. II, we shall hence briefly define the model that is studied and the quantities that will be analyzed and comment on the simulation methods. Sec. III briefly reviews the scaling predictions, including the finite-size scaling results for the case where both $`D`$ and the lateral linear dimension $`L`$ are finite. Sec. IV then presents our results on $`T_c(D)`$ and the critical fields $`H_c(D)`$, on which our tests of Eqs. (1), (2) are based. Sec. V discusses those aspects of our results which are pertinent to a test of two-dimensional criticality, Eq. (3), while Sec. VI summarizes our conclusions.

## II Model and simulation technique

Invoking the well-known isomorphism between the lattice gas model of fluids and the Ising model of magnetism (see e.g. Ref.  for details), we study the Ising model on the simple cubic lattice in the presence of a bulk field $`H`$ and a surface field $`H_1`$,

$$=J\underset{i,j}{}S_iS_jH\underset{i}{}S_iH_1\underset{iϵ\mathrm{surfaces}}{}S_i,S_i=\pm 1\text{,}$$
(4)
where the exchange interaction $`J`$ is only present between nearest neighbors on the lattice. Note that phase coexistence in the bulk (phases with positive and negative magnetization correspond to gas and liquid phases of the fluid, respectively) corresponds to $`H=0`$, see Fig. 1. In the thin film, one trivially obtains the result that for zero temperature phase coexistence occurs for  $`H_{\text{coex}}(D,T=0)=2H_1/D`$, but for $`T>0`$ the variation of $`H_{\text{coex}}(D,T)`$ is nontrivial. While previous work  was mostly interested in the behavior of $`H_{\text{coex}}(D,T)`$ near the temperature $`T_w(H_1)`$ of the wetting transition , we consider here scaled surface fields $`H_1D^{\mathrm{\Delta }_1/\nu }`$ small enough such that we stay in the nonwet regime of the surface phase diagram of the semi-infinite system  throughout, although we consider the vicinity of $`T_c^{\mathrm{}}=T_c(D=\mathrm{})`$. Measuring all lengths in units of the lattice spacing, we consider film thicknesses $`D=4,8,12,16,24,28`$ and $`32`$ for an $`L\times L\times D`$ geometry, varying $`L`$ over an as wide range as is practical, from the point of view of available computer ressources. In the $`x,y`$-directions parallel to the thin film, we apply periodic boundary conditions as usual .

In order to be able to find $`T_c(D)`$ and $`H_c(D)H_{\text{coex}}(D,T=T_c(D))`$ reliably, we have to use aspect ratios $`L/D1`$. Although the choices of film thickness as quoted above are not extremely large, it is clear that use of rather large linear dimensions $`L`$ is mandatory for obtaining reliable results. If we would use the Metropolis algorithm , as done in Refs. , or the heatbath algorithm , “critical slowing down”  would be a serious problem: i.e., the “time” $`\tau `$ over which subsequently generated system configurations are correlated varies as

$$\tau L^z\text{ with }z2.16(d=2)\text{ or }z2.03(d=3),$$
(5)
the prefactor in this power law being of order unity if Monte Carlo time is measured in units of attempted Monte Carlo steps (MCS) per spin. Since we wish to use linear dimensions of the order of $`L10^2`$, relaxation times of the order of $`10^4`$ MCS easily result. Given the fact that quantities like the specific heat $`C_v`$ and the susceptibility $`\chi `$, recorded from fluctuations of energy and magnetization,

$$C_v=\left(^2^2\right)/(L^2Dk_\text{B}T^2),\chi =\left\{\left(\underset{i}{}S_i\right)^2\underset{i}{}S_i^2\right\}/(L^2Dk_\text{B}T),$$
(6)
are non-self-averaging , one needs $`n1`$ statistically independent observations (i.e., separated by time intervals $`\mathrm{\Delta }t>\tau `$) to obtain $`C_v`$ and $`\chi `$ with small enough error (the relative error of these quantities is  $`\sqrt{2/n}`$, irrespective of $`L`$ and $`D`$). For this reason, it is clear that the use of cluster algorithms which reduce critical slowing down  is highly desirable. However, for the present problem where both a bulk magnetic field and a surface magnetic field of competing sign are present {Eq. (4)} application of cluster algorithms is nontrivial. It turns out that an extension of the “ghost spin algorithm”  to the present problem is rather straightforwardly possible . The coupling of spins to a magnetic field is thereby treated as if it were an additional infinite-range exchange coupling to an extra spin $`S_G=\pm 1`$. This coupling has the strength $`h=|H|`$ for spins in the interior of the film and $`h=|H_1+H|`$ for spins in the surface layers. In addition to putting bonds in clusters of spins (inside a cluster all spins are connected by exchange interactions and have the same sign) with probability  $`p_B=1\mathrm{exp}(2J/k_\text{B}T)`$ one also puts bonds between the spins in clusters and the ghost spin $`p_G=1\mathrm{exp}(2h/k_\text{B}T)`$, if the orientation of the spins in the cluster is the same as that of the ghost spin {which is $`S_G=`$sign $`(H)`$ for interior spins and $`S_G=`$sign $`(H_1+H)`$ for spins in the surface planes, respectively}.

While this extension of the cluster algorithm to the case of nonzero bulk and surface fields is formally exact, discussion of its efficiency is a rather delicate problem: in fact, if $`h/k_\text{B}T`$ is of order unity, also $`p_G`$ is of order unity and the infinite-range character of this coupling then implies that huge clusters containing a large fraction of the entire simulation volume would be created most of the time! It is clear that under such circumstances the algorithm would be very inefficient; as in the case of zero field it is necessary for a good performance of a cluster algorithm that typically large clusters are created but a single large cluster must contain only a negligible fraction of the total volume in the thermodynamic limit. As a consequence, one needs $`h/k_\text{B}T1`$, and since $`k_\text{B}T`$ is in the range of $`3.54.5`$ we have thus chosen to work with a single value of the surface field, namely $`H_1=0.015J`$. Even for this small value – note that the corresponding value of $`H`$ is typically one or two orders of magnitude smaller, see below – the performance of the algorithm has significantly deteriorated, in comparison with the case without any magnetic fields. This fact can be clearly demonstrated by a binning analysis  of the magnetization $`m`$ in the system: the $`N`$ (dynamically correlated) subsequent observations $`m_\nu =(1/N)\underset{i=1}{\overset{N}{}}S_i^\nu `$ are grouped into $`n=N/N_b`$ blocks of length $`N_b`$, from which block averages $`\stackrel{~}{m}_\mu `$ of the corresponding $`N_b`$ observations {$`m_\nu `$} belonging to the block with index $`\mu `$ are formed. Then

$$\mathrm{\Delta }m[n(n1)]^{1/2}\sqrt{\underset{i=1}{\overset{n}{}}(\stackrel{~}{m}_\mu \overline{m})^2}$$
(7)
(where $`\overline{m}=n^1\underset{i=1}{\overset{n}{}}\stackrel{~}{m}_i)`$ is studied as a function of block length $`N_b`$ (Fig. 2): when $`\mathrm{\Delta }m`$ is independent of $`N_b`$, the subsequent $`\stackrel{~}{m}_\mu `$ are statistically independent, and $`\mathrm{\Delta }m`$ is a good estimation of the statistical error; otherwise one sees a systematic increase of $`\mathrm{\Delta }m`$ with $`N_b`$ and the value of $`N_b`$ needed to reach a saturation value yields an estimate for the correlation time. For the Metropolis algorithm and the chosen system size ($`L`$=128), even for $`N_b=4000`$ one is far from saturation, and hence it is clear that this algorithm would be very impractical for the present problem. For the cluster algorithm and $`H_1=0`$, on the other hand, $`\mathrm{\Delta }m`$ vs. $`N_b`$ is essentially constant, $`\mathrm{\Delta }m0.0004`$, the correlation time being of order unity, as expected . However, this is not so for the cluster algorithm in the case $`H_1=0.015`$: $`\mathrm{\Delta }m`$ saturates at a plateau of about $`\mathrm{\Delta }m0.007`$, i.e. the error is almost a factor 20 larger, and the correlation time is of the order of $`\tau _m280`$ Monte Carlo steps in the example shown in Fig. 2(b). Thus, while the gain of the cluster algorithm in the zero field case compared to the Metropolis algorithm is very significant, in our problem it is only rather modest! This - somewhat unexpected - dramatic decrease of the efficiency of the cluster algorithm with increasing strength of the surface (and bulk) fields has prevented us both from studying systems larger than $`L=128`$ and from studying the dependence on $`H_1`$ systematically. Runs of length up to 1.2 million Monte Carlo steps (MCS) were performed.

As is well known , cluster algorithms at critical points of Ising systems are rather sensitive to correlations among the pseudorandom numbers produced by the random number generator. In the present work, we used thus an improved version of the standard ”R250” generator , where two versions {one based on the pair of integers (250, 103) the other with the pair (521, 168)} are combined with the logical XOR operation.

In order to make best use of our simulation data, we apply standard multihistogram interpolation techniques . Note that we needed a three-dimensional histogram $`P(E,m,m_1)`$, $`E`$ being the exchange energy, $`m_1`$ the magnetization in the surface plane, in order to allow reweightings in the full parameter space of independent control variables ($`T,H,H_1`$) and hence the storage requirements for $`P`$ are nontrivial. However, noting that all measurements of $`E,m,m_1`$ can be represented by integers, each integer needing 4 Byte, we can store the time series of 106 observations with a storage of 12 MByte, irrespective of the choices of $`L`$ and $`D`$.

The multihistogram reweighting with respect to the bulk field $`H`$ is crucial in order to be able to find the field $`H_{\text{coex}}(T)`$, along which for $`T<T_c(D)`$ two-phase coexistence occurs, applying the “equal weight rule” : In the space of variables $`(E,m,m_1)`$, the two phases show up as separate peaks of $`P(E,m,m_1)`$ {or $`P(m,m_1)`$, respectively, see Fig. 3(a), when one studies an isotherm one can integrate out $`E`$, of course}, which have precisely the same weight at $`H=H_{\text{coex}}(T)`$ while for $`HH_{\text{coex}}(T)`$ (but not too far away from it) the two peaks can still be identified but have different weight. With the multihistogram reweighting, a small number of simulation points suffices to generate the curve $`H_{\text{coex}}(T)`$ {and its extrapolation into the regime $`T>T_c(D)`$} with reasonable precision, see Fig. 3(b). Since near $`T_c(D)`$ the free energy differences between the two phases are very small also off coexistence, the statistical error in the estimation of $`H_{\text{coex}}(T)`$ is not negligible, and also systematic errors, since $`L/D`$ is not large enough, need to be considered. The latter problem also affects the estimation of $`T_c(D)`$, as will be discussed in Sec. IV.

## III SCALING PREDICTIONS

For completeness, we first summarize the pertinent predictions of the scaling theory for thin Ising films near the critical point , assuming the lateral linear dimension $`L`$ infinite, and consider the extension  to finite $`L`$ in the following. The singular part of the free energy per spin is assumed to scale as follows

$$f_{\text{sing}}(D,T,H,H_1)|t|^{2\alpha }\stackrel{~}{f}_\pm (D|t|^\nu ,H|t|^{\mathrm{}},H_1|t|^\mathrm{}_1),$$
(8)
where $`\alpha `$ is the exponent of the specific heat of the three-dimensional Ising model, $`t=\left(TT_c\left(\mathrm{}\right)\right)/T_c\left(\mathrm{}\right),\stackrel{~}{f}_\pm `$ is a “scaling function” (with two different branches, referring to the sign of $`t`$), and the other exponents have already been defined in Sec. I.

Now it is convenient to introduce the scaling variables

$$xD|t|^\nu ,wH_1D^{\mathrm{\Delta }_1/\nu },$$
(9)
and then Eq. (8) can also be written as, eliminating $`|t|`$ from the arguments of $`\stackrel{~}{f}_\pm `$,

$$f_{\text{sing}}(D,T,H,H_1)|t|^{2\alpha }\stackrel{~}{f}_\pm (x,\frac{HD^{\mathrm{\Delta }/\nu }}{x^{\mathrm{\Delta }/\nu }},\frac{w}{x^{\mathrm{\Delta }_1/\nu }}).$$
(10)
Since the critical point of the thin film is shifted relative to the bulk critical point $`T_c\left(\mathrm{}\right)`$, it must correspond to a singular behavior of the scaling function $`\stackrel{~}{f}_\pm `$. At fixed $`H_1`$ and fixed $`D`$ this means the scaling function $`\stackrel{~}{f}_\pm `$ $`\left(x,y,y_1=w/x^{\mathrm{\Delta }_1/\nu }\right)`$ has a singularity at a point $`x_c\left(w\right),`$ $`y_c\left(w\right)`$. Therefore the shifts $`\mathrm{\Delta }T_c\left(D\right),\mathrm{\Delta }H_c\left(D\right)`$ follow as

$$\mathrm{\Delta }T_c=T_c(D,H_1)T_c\left(\mathrm{}\right)=B_TD^{1/\nu }X_c\left(CH_1D^{\mathrm{\Delta }_1/\nu }\right),$$
(11)
$$\mathrm{\Delta }H_cH_c(D,H_1)=B_HD^{\mathrm{\Delta }/\nu }Y_c(CH_1D^{\mathrm{\Delta }_1/\nu }).$$
(12)
The scaling functions $`X_c,Y_c`$ are universal, while $`B_T,B_H`$ and $`C`$ are non-universal critical amplitudes, which are normalized such that $`X_c(Cw)=1+(Cw)^2+\mathrm{},`$ $`Y_c(Cw)Cw+o(Cw)^3`$. Note that both functions are analytic for $`w0`$, and $`\mathrm{\Delta }T_c`$ should be an even function of $`H_1`$ and hence $`w`$, while $`\mathrm{\Delta }H_c`$ must be an odd function of $`H_1`$. From these considerations, for small $`H_1`$ Eqs. (1), (2) result, remembering  that $`\mu _c(D)\mu _c(\mathrm{})=2H_c`$.

An alternative argument for Eq. (1), which also elucidates how this equation fits together with the Kelvin equation $`(H_cH_1/D)`$ in the critical region for large enough $`D`$, derives from a consideration of phase coexistence for temperatures slightly below $`T_c(D)`$. If $`H_1=H=0`$, we would have two coexisting phases with magnetization profiles $`m^+(z)`$ and $`m^{}(z)=m^+(z)`$ across the film, and both states have the same free energy $`F_+(0,0)=F_{}(0,0)`$. Since these profiles are smooth functions of $`H`$ and $`H_1`$, an expansion of the free energies around $`F_+(0,0)`$ {or $`F_{}(0,0)`$, respectively} yields

$$\mathrm{\Delta }F_+F_+(0,0)F_+(H,H_1)=\overline{m^+}HDL^2+2m_1^+H_1L^2,$$
(13)
$$\mathrm{\Delta }F_{}F_{}(0,0)F_{}(H,H_1)=\overline{m^{}}HDL^2+2m_1^{}H_1L^2,$$
(14)
where $`\overline{m^+},\overline{m^{}}`$ refer to the average over the magnetization profile in the respective states, and $`m_1^+,m_1^{}`$ the layer magnetizations in the surface layer. To leading order for small $`H_1`$ and small $`H`$ in Eqs. (13,14), $`\overline{m^+},\overline{m^{}}`$ and $`m_1^+,m_1^{}`$ are to be taken at zero fields, and thus satisfy the symmetry $`\overline{m^{}}=\overline{m^+}`$ and $`m_1^{}=m_1^+`$. Phase coexistence occurs for $`\mathrm{\Delta }F^+=\mathrm{\Delta }F^{}`$, and hence the Kelvin equation results,

$$H_{\text{coex}}(D,T,H_1)\frac{2H_1}{D}\frac{m_1^+(D,T)}{\overline{m^+}(D,T)}.$$
(15)
Assuming then that $`D\xi `$, the correlation length in the bulk, $`m^+(z)`$ will approach the bulk spontaneous magnetization $`m_b=\widehat{B}_b(t)^\beta `$ almost everywhere, and hence $`\overline{m_+}\widehat{B}_b(t)^\beta `$, $`\widehat{B}_b`$ being the respective critical amplitude. Likewise $`m_1^+(D,T)`$ approaches the surface layer magnetization of a semi-infinite system  $`m_1=\widehat{B}_1(t)^{\beta _1}`$, with $`\widehat{B}_1`$ the corresponding amplitude. Therefore Eq. (15) becomes in this limit

$$H_{\text{coex}}(D,T,H_1)\frac{2H_1}{D}\frac{\widehat{B}_1}{\widehat{B}_b}(t)^{\beta _1\beta },x=D|t|^\nu \mathrm{}.$$
(16)
Since  $`\beta 0.325`$ and  $`\beta _10.780.80`$, we see that the coefficient of the term $`H_1/D`$ in Eq. (16) gets smaller and smaller as $`|t|`$ gets smaller. Due to this vanishing coefficient in the limit $`|t|0`$ a smooth crossover between the Kelvin equation, Eqs. (15, 16), and Eq. (1) becomes possible. Remembering that for $`D`$ finite there is a shift of $`T_c`$ as given by Eq. (2), we can further conclude that the critical field should be of the order

$$H_c(D_1,H_1)H_{\text{coex}}(D,T_c(D),H_1)\frac{2H_1}{D}D^{(\beta _1\beta )/\nu }=2H_1D^{(\mathrm{\Delta }\mathrm{\Delta }_1)/\nu },$$
(17)
where in the last step the standard scaling relations $`\beta =2\alpha \mathrm{\Delta },\beta _1=2\alpha \nu \mathrm{\Delta }_1`$ were used. Eq. (17) obviously is nicely compatible with Eq. (1).

Eq. (16) allows an interesting conclusion to be drawn on the slope of the coexistence curve at the critical point of the film (cf. Fig. 1). We find first for the angle $`\theta (t)`$ describing this slope for $`T<T_c(D)`$

$$\mathrm{tan}(\theta )(H/T)_{H_1}=\frac{2H_1}{T_c(\mathrm{})}\frac{\beta _1\beta }{D}\frac{\widehat{B}_1}{\widehat{B}_b}(t)^{\beta _1\beta 1},x=D|t|^\nu \mathrm{}.$$
(18)
Since $`\beta _1<1+\beta `$, the exponent of $`(t)`$ is negative, and thus for the considered limit the slope diverges (i.e., varying $`t`$ at very large but fixed $`D`$). However, this limit does not include the limiting slope at $`TT_c(D)`$ itself, since then Eqs. (2), (11) yield $`tD^{1/\nu }`$, and hence

$$\mathrm{tan}(\theta )H_1D^{(\beta _1\beta 1)/\nu 1}=H_1D^{(\mathrm{\Delta }\mathrm{\Delta }_11)/\nu }.$$
(19)
In Landau theory, $`\mathrm{\Delta }=3/2,\mathrm{\Delta }_1=1/2`$ and hence the power of $`D`$ vanishes, i.e., the slope is nonzero and finite at $`T_c(D)`$ in the limit $`D\mathrm{}`$. For the three-dimensional Ising model, the best available exponent estimates  imply $`(\mathrm{\Delta }\mathrm{\Delta }_11)/\nu 0.120.16`$, i.e. $`\theta 0`$ for $`D\mathrm{}`$! This result also implies that for capillary condensation field mixing effects  are asymptotically not very important.

We now briefly consider the extension of the scaling theory to include finite-size effects due to the finite lateral linear dimension $`L`$ {in Eqs. (13, 14) we have assumed the limit $`L\mathrm{}`$ throughout}. This can simply be done by including the aspect ratio $`L/D`$ as an additional scaling variable in Eqs. (8, 9), which we then rewrite as follows

$$f_{\text{sing}}(D,T,H,H_1,L)D^3\stackrel{~}{f}(D^{1/\nu }t,L/D,HD^{\mathrm{\Delta }/\nu },H_1D^{\mathrm{\Delta }_1/\nu }).$$
(20)
Since for finite $`L`$ the free energy and its derivatives are smooth functions of $`t`$, it is more convenient to use $`D^{1/\nu }t`$ rather than $`x=D|t|^\nu `$ as a scaling variable. From Eq. (20), we immediately obtain the following scaling results for the specific heat, the magnetization and the susceptibility of the thin film

$$C_v=D^{\alpha /\nu }\stackrel{~}{C}(D^{1/\nu }t,L/D,HD^{\mathrm{\Delta }/\nu },H_1D^{\mathrm{\Delta }_1/\nu }),$$
(21)
$$m=D^{\beta /\nu }\stackrel{~}{m}(D^{1/\nu }t,L/D,HD^{\mathrm{\Delta }/\nu },H_1D^{\mathrm{\Delta }_1/\nu }),$$
(22)
$$\chi =D^{\gamma /\nu }\stackrel{~}{\chi }(D^{1/\nu }t,L/D,HD^{\mathrm{\Delta }/\nu },H_1D^{\mathrm{\Delta }_1/\nu }),$$
(23)
where $`\stackrel{~}{C},\stackrel{~}{m}`$ and $`\stackrel{~}{\chi }`$ are appropriate scaling functions. Since we choose $`H_1`$ fixed, $`D`$ fixed, $`H_1D^{^{\mathrm{\Delta }_1/\nu }}1`$, and $`H`$ is chosen according to Eqs. (16, 17) {in practice this is done by applying the reweighting technique and the equal area rule, cf. Fig. 3} the last two arguments $`HD^{\mathrm{\Delta }/\nu },H_1D^{\mathrm{\Delta }_1/\nu }`$ in Eqs. (21)–(23) can be ignored in the following discussion.

Now in the limit $`L\mathrm{}`$ we expect that $`C_v`$ exhibits a logarithmic singularity for $`TT_c(D)`$, while the critical part of the magnetization $`m_{\text{crit}}mm(T_c(D),H,H_1)`$ should behave as

$$m_{\text{crit}}[1T/T_c(D)]^{\beta _2},\beta _2=1/8,$$
(24)
and the susceptibility

$$\chi |1T/T_c(D)|^{\gamma _2},\gamma _2=7/4.$$
(25)
For finite $`L`$, however, these singularities are all rounded off and we rather expect that both $`C_v`$ and $`\chi `$ exhibit maxima of finite height at temperatures $`T_{\mathrm{max}}^c`$ $`\left(D\right)`$, $`T_{\mathrm{max}}^\chi \left(D\right)`$. From Eqs. (21), (23), we readily predict (in the limit $`|H_1|D^{\mathrm{\Delta }_1/\nu }1`$)

$$\frac{T_{\mathrm{max}}^c\left(D\right)T_c(D)}{T_c(\mathrm{})}D^{1/\nu }=\mathrm{\Delta }\stackrel{~}{T}_{\mathrm{max}}^c\left(L/D\right),$$
(26)
$$\frac{T_{\mathrm{max}}^\chi \left(D\right)T_c(D)}{T_c(\mathrm{})}D^{1/\nu }=\mathrm{\Delta }\stackrel{~}{T}_{\mathrm{max}}^\chi \left(L/D\right),$$
(27)
with $`\mathrm{\Delta }\stackrel{~}{T}_{\mathrm{max}}^c\left(L/D\right)`$, $`\mathrm{\Delta }\stackrel{~}{T}_{\mathrm{max}}^\chi \left(L/D\right)`$ being suitable scaling functions that describe the shift of these maxima positions as functions of the aspect ratio $`L/D`$. From this analysis one also can predict  how the height of the maxima should depend on $`D`$ and $`L`$, for $`L>>D`$

$$C_v^{\mathrm{max}}D^{\alpha /\nu }\mathrm{ln}(L/D),$$
(28)
$$\chi ^{\mathrm{max}}D^{\gamma /\nu 7/4}\text{ }L^{7/4}\text{,}$$
(29)
and how the absolute value of the order paramter should decrease at $`T_c(D)`$,

$$|m_{\text{crit}}|_{T=T_c\left(D\right)}D^{1/8\beta /\gamma }L^{1/8}.$$
(30)
Finally, in the limit $`L\mathrm{}`$ the $`D`$-dependence of the critical amplitudes associated with the two-dimensional critical behavior (see Ref.  for a more detailed discussion), defining now $`\stackrel{~}{t}=\left[TT_c\left(D\right)\right]/T_c\left(\mathrm{}\right)`$, can be read off from the following equations,

$$C_vD^{\alpha /\nu }\mathrm{ln}|\stackrel{~}{t}|,$$
(31)
$$m_{\text{crit}}D^{\left(1/8\beta \right)/\nu }\left(\stackrel{~}{t}\right)^{1/8},$$
(32)
$$\chi D^{\left(\gamma 7/4\right)/\nu }|\stackrel{~}{t}|^{7/4}\text{ }.$$
(33)
Due to the crossover scaling between two- and three-dimensional critical behavior, that Eqs. (20), (21), (22), and (23) describe, a singular dependence of the various critical amplitudes on film thickness results at the capillary condensation critical point.

## IV Numerical Results on $`T_c\left(D\right)`$ and $`H_c\left(D\right).`$

For locating critical points in the bulk, a convenient method is to study the fourth order cumulant of the order parameter

$$U_L\left(T\right)=1m_{\text{crit}}^4/\left[3m_{\text{crit}}^2^2\right]$$
(34)
for a range of linear dimensions $`L`$ as a function of temperature, and to look for a common intersection point  {which for the universality class of the two-dimensional Ising model, should have the value  $`U^{}=U_L\left(T=T_c\right)`$ $`=0.610690\left(1\right)`$}. In our case, we have to follow a path along $`H=H_{\text{coex}}\left(T\right)`$ {as shown in Fig. 3(b)} when we record these cumulants, and since this path is not exactly known but only within some numerical error, it is clear that this method is more difficult to apply than for ordinary bulk Ising models. In addition, even for small $`D`$ the data are plagued by crossover scaling effects (Fig. 4): the curves for the values of $`L`$ that are practically available do not intersect in a common point, but rather the intersection points are scattered and fall below the theoretical value $`U^{}`$. This failure of verifying the common intersection points is not unexpected, since Fig. 4 includes data for which the aspect ratio $`L/D`$ is as small as 4 (a) or even 2 (b), rather than only data for which $`L/D1`$. In fact, from the treatment of the previous section we can readily conclude that

$$U_L\left(T=T_c\left(D\right)\right)=\stackrel{~}{U}\left(L/D\right)$$
(35)
and only in the limit $`L/D\mathrm{}`$ shall we have $`U\left(\mathrm{}\right)=U^{}`$.

An alternative and widely used recipe to find the critical temperature is to try an extrapolation of the maxima of the specific heat and susceptibility versus $`L^{1/\nu }`$or of the cumulant intersection points. Considering the intersection of $`U_L\left(T\right)`$ and $`U_{bL}\left(t\right)`$ with a scale factor $`b>1`$, it can be argued  that corrections to finite-size scaling lead to a shift of the intersection point that varies with $`b`$ proportional to $`\left[b^{1/\nu }1\right]^1`$ for large $`b`$. Fig. 5 shows some attempts to carry out such extrapolations, again for $`D=8`$ and $`D=28`$ (data for all other choices of $`D`$ can be found in  and look similar). These figures show that $`T_{\mathrm{max}}^c`$ approaches $`T_c\left(D\right)`$ in a non-monotonic fashion, and also the curve $`T_{\mathrm{max}}^\chi `$ vs. $`L^{1/\nu }`$ is distinctly non-linear. Fitting asymptotic straight lines to both data sets one obtains results for $`T_c\left(D\right)`$ that are roughly compatible with each other, and with the (linear) extrapolation of the cumulant intersections. Although the accuracy of $`T_c\left(D\right)`$ obtained in this way is several orders of magnitude less than in the case of the bulk three-dimensional Ising model , the data are accurate enough to allow a useful test of Eqs. (1), (2).

The consistency of our analysis can be checked further by testing for the scaling behavior predicted in Eqs. (26), (27), see Fig. 6. Here all data points are included for all values of $`D`$ and $`L`$ that have been studied and $`T_c(D)`$ is chosen such that the best data collapse is achieved. It is seen that the non-monotonic variation of the temperature at which the specific heat has its maximum is an intrinsic property of this scaling function describing the system shape effects in terms of the aspect ratio $`D/L`$ of the simulation box. The interpolating curves are simple parabolic fits which translate back into the solid lines in the left part of Fig. 5.

The values of $`T_c\left(D\right)`$ that we have determined as shown in Figs. 5, 6 are collected in Table I, which includes also our estimates for $`H_c\left(D\right)`$. Log-log plots of these data versus $`D`$ almost look like straight lines, however, there is a slight but systematic curvature, and if this curvature were disregarded and straight lines were fitted to all the data nevertheless, the resulting effective exponents would systematically deviate from the theoretial predictions in Eqs. (1), (2).

Better results are obtained if one fits effective exponents from successive thicknesses only $`\left(D=4,8,12;D=8,12,16;\mathrm{};D=24,28,32\right)`$, which can be extrapolated vs. $`1/D`$ reasonably well, and converge nicely towards the theoretical predictions (Fig. 7), namely $`1/\nu 1.587`$ and $`\left(\mathrm{\Delta }\mathrm{\Delta }_1\right)/\nu 1.75`$. Conversely, if Fig. 7(b) was taken as an independent estimation of the exponent $`\mathrm{\Delta }_1`$, we would obtain $`\mathrm{\Delta }_1=0.459(13)`$, which indeed is compatible with the existing recent estimates .

## V A Test of Two-Dimensional Critical Behavior

In this section, we are concerned with the question to what extent the data provide some evidence for the prediction (Eq. 3) that the capillary condensation critical point displays critical exponents of the two-dimensional Ising universality class. Since the accessible values of the lateral linear dimension $`L`$ are not very large, however, we cannot expect that a regime can be reached where the parallel correlation length $`\xi _{||}`$ satisfies the criterion $`D\xi _{||}L`$ \- only in such a regime a direct observation of these power laws would be possible. Hence we attempt to study the critical behavior again via a finite-size scaling analysis, using ,

$$|m_{\text{crit}}|L^v=\stackrel{~}{M}\left(L^u\stackrel{~}{t}\right),$$
(36)
$$\chi L^w=\stackrel{~}{\chi }\left(L^u\stackrel{~}{t}\right),$$
(37)
where the exponents $`u,v,w`$ should take the values

$$u=1/\nu _2=1,v=\beta /\nu _2=1/8,w=\gamma _2/\nu _2=7/4.$$
(38)
Eqs. (36), (37), and (38) are appropriate if $`D\xi _{||}`$ still holds but $`\xi _{||}`$ and $`L`$ are of the same order. In practice, however, also the condition $`D\xi _{||}`$ is hard to satisfy since we wish to include some data for which $`L/D`$ is not very large. It then helps to relax the theoretical condition, Eq. (38), and rather treat $`u,v,w`$ as effective exponents : in this way, one can take into account to some extent the corrections to finite-size scaling arising from the crossover between two- and three-dimensional Ising critical behavior.

Fig. 8 shows that this procedure works reasonably well, and Table II gives a listing of the fit exponents $`u,v,w`$, and corresponding effective exponents $`\nu _{\text{eff}}=1/u`$, $`\beta _{\text{eff}}=v/u`$ and $`\gamma _{\text{eff}}=w/u`$. It is seen from Table II that both $`u,v`$ and $`w`$ gradually change from the two-dimensional values towards the three-dimensional ones, although even for $`D=32`$ one is still far away from the theoretical values for the latter. While $`\beta _{\text{eff}}`$ has increased significantly, $`\gamma _{\text{eff}}`$ within the accuracy of this estimation has hardly changed at all. If we consider an effective dimensionality from the hyperscaling relation , defined as $`d_{\text{eff}}=(\gamma _{\text{eff}}+2\beta _{\text{eff}})/\gamma _{\text{eff}}=w+2v`$, we find $`d_{\text{eff}}=2.0\pm 0.15`$, and there is no systematic trend with $`D`$. While the latter observation is in accord with a previous study using “neutral walls” , where $`H_c(D)0`$, we have obtained in the present work a much better evidence that for small $`D`$ the behavior is compatible with two-dimensional Ising criticality than was possible in the latter model . Note also that in the present study there is a rather broad range of $`D`$ where $`\nu _{\text{eff}}>1`$, which was not the case in . Due to the systematic problems of fitting several effective exponents from somewhat noisy data and the restricted range over which this fit is applicable we do not think that these discrepancies are a proof of non-universal crossover behavior, however. We feel that this problem needs a more careful study.

## VI CONCLUSIONS

In this paper Monte Carlo simulations have been presented attempting to test theoretical predictions about the capillary condensation critical point. Using an extension of the Swendsen-Wang cluster algorithm including competing surface and bulk magnetic fields, for Ising films with thicknesses $`D=4,8,12,16,24,28`$ and $`32`$ the critical temperature $`T_c(D)`$ and corresponding critical field $`H_c(D)`$ for a surface magnetic field $`H_1`$ have been estimated. The data are compatible with the power laws presented about 20 years ago by Fisher and Nakanishi. Also the expected two-dimensional critical behavior is compatible with our data, though the accuracy of the resulting effective exponents is rather low (Table II) and hence a more convincing proof would be desirable, but is not feasible with the present computer ressources.

A challenging problem that we have not solved is the development of an efficient version of the cluster algorithm that allows to work with surface and bulk fields that are not extremely small. The algorithm that we have used was much less efficient even for $`H_1=0.015J`$ than for $`H_1=0`$, and a study of capillary condensation critical points over the range where $`(H_1/J)D^{\mathrm{\Delta }_1/\nu }`$ is not small, and hence the nonlinear part of the scaling function $`Y_c(CH_1D^{\mathrm{\Delta }_1/\nu })`$ would be probed, turned out not to be feasible either. Thus, in spite of a longstanding effort to deal with theory and simulation of capillary condensation there remain still some missing links. A particularly intriguing problem is to elucidate the crossover between three-dimensional and two-dimensional critical behavior in these thin films. Finally, it is hoped that the present study provides an incentive to address this problem also by suitable experiments.

Acknowledgements: One of us (O. D.) acknowledges support from the Deutsche Forschungsgemeinschaft (DFG) under grant No. Bi314/16, another (W. J.) acknowledges support from the Heisenberg program of the DFG and from Sonderforschungsbereich 262/D1.

TABLE I: Critical temperatures and fields

| $`D`$ | $`T_c(D)`$ | $`H_c(D)/J`$ |
| --- | --- | --- |
| 4 | 3.8705(3) | 0.006644(32) |
| 8 | 4.2409(3) | 0.002528(14) |
| 12 | 4.3561(3) | 0.001367(10) |
| 16 | 4.4084(3) | 0.000867(6) |
| 24 | 4.4549(5) | 0.000448(3) |
| 28 | 4.4665(4) | 0.000348(3) |
| 32 | 4.4749(5) | 0.000279(3) |

| $`\mathrm{}`$ | 4.51152(2) | 0 |
| --- | --- | --- |

TABLE II: Effective exponents for order parameter and susceptibility

| $`D`$ | $`u`$ | $`v`$ | $`w`$ | $`\beta _{\text{eff}}`$ | $`\gamma _{\text{eff}}`$ | $`\nu _{\text{eff}}`$ |
| --- | --- | --- | --- | --- | --- | --- |
| 2-dim | 1 | 1/8 | 7/4 | 1/8 | 7/4 | 1 |
| 4 | 0.956 | 0.126 | 1.67 | 0.132 | 1.75 | 1.064 |
| 8 | 1.018 | 0.136 | 1.72 | 0.133 | 1.69 | 0.982 |
| 12 | 0.944 | 0.138 | 1.67 | 0.146 | 1.77 | 1.059 |
| 16 | 0.938 | 0.139 | 1.61 | 0.148 | 1.72 | 1.066 |
| 24 | 0.898 | 0.145 | 1.53 | 0.161 | 1.70 | 1.114 |
| 28 | 0.853 | 0.141 | 1.48 | 0.165 | 1.74 | 1.172 |
| 32 | 0.884 | 0.155 | 1.54 | 0.175 | 1.74 | 1.131 |
| 3-dim | 1.587 | 0.518 | 1.96 | 0.327 | 1.24 | 0.630 |


format: text

text:
# Phenomenology of neutrino physics in the Kaluza-Klein theories of low scale gravity

## Abstract

We discuss the phenomenological consequences of theories which describe sterile neutrinos in large extra dimensions. We show that the Kaluza-Klein tower of the singlet neutrinos, albeit tiny individual contribution in electroweak processes, act cumulatively, giving rise to non-universality of the weak interactions of the light neutrinos and to flavour-violating radiative processes. Owing to these non-decoupling effects of th Kaluza–Klein neutrinos, we derive strong constraints on the parameters of the theory that originates from the non-observation of flavour-violating and universality-breaking phenomena. In this theory we propose a four-neutrino model which can reconcile the existing data coming from underground experiments in terms of neutrino oscillations, together with the hint from the LSND experiment and a possible neutrino contribution to the hot dark matter of the Universe.

Recently, it has been proposed  the radical possibility that the observed small value of the gravitational constat at long distance is ascribed to the spreading of the gravitational force in $`\delta `$ large extra spatial dimensions. It follows that the Standard Model (SM) fields are confined to a 3-brane configuration, while the large compactified dimensions are probed only by gravity and singlet, under SM gauge group, fields.

In this presentation we examine the phenomenological consequences of higher - dimensional isosinglet neutrinos on collider and lower energy experiments.

For our phenomenological study we adopt a variant  of the model discussed in Ref.  . For definiteness, we consider a model that minimally extends the SM-field content by one singlet Dirac neutrino, $`N(x,y)`$, which propagates in a $`[1+(3+\delta )]`$-dimensional Minkowski space. The $`y`$-coordinates are compactified on a circle of radius $`R`$ by applying the periodic identification: $`yy+2\pi R`$.

The relevant part of the Lagrangian of our minimal model reads

$`{\displaystyle \underset{0}{\overset{2\pi R}{}}}d\stackrel{}{y}[\overline{N}(i\gamma ^\mu _\mu +i\gamma _\stackrel{}{y}_\stackrel{}{y}m)N`$
$`+\delta (\stackrel{}{y})({\displaystyle \underset{l=e,\mu ,\tau }{}}{\displaystyle \frac{h_l}{M_F^{\delta /2}}}L_l\stackrel{~}{\mathrm{\Phi }}\xi +\mathrm{H}.\mathrm{c}.)],`$ (1)
where $`\stackrel{~}{\mathrm{\Phi }}`$ ( $`<\stackrel{~}{\mathrm{\Phi }}>=v=174`$Gev ) and $`L_l`$ are higgs and lepton doublets, respectively. We have assumed that only one two-component spinor $`\xi `$ from the high dimensional spinor $`N`$ couples to our 3 brane.

We can now express the two-component spinor $`\xi `$ and its dirac partner $`\eta `$ from high dimensional spinor $`N`$ in terms of a Fourier series expansion as follows:

$`[\xi (x,y),\eta (x,y)]=`$ (2)
$`{\displaystyle \frac{1}{(2\pi R)^{\delta /2}}}{\displaystyle \underset{\stackrel{}{n}}{}}[\xi _\stackrel{}{n}(x),\eta _\stackrel{}{n}(x)]\mathrm{exp}\left({\displaystyle \frac{i\stackrel{}{n}\stackrel{}{y}}{R}}\right).`$
Substituting into the Eq. (Phenomenology of neutrino physics in the Kaluza-Klein theories of low scale gravity) and then performing the $`\stackrel{}{y}`$ integration yields

$`{\displaystyle \underset{\stackrel{}{n}}{}}\{\overline{\xi }_\stackrel{}{n}(i\overline{\sigma }^\mu _\mu )\xi _\stackrel{}{n}+\overline{\eta }_\stackrel{}{n}(i\overline{\sigma }^\mu _\mu )\eta _\stackrel{}{n}`$ (3)
$`[\xi _\stackrel{}{n}((m+{\displaystyle \frac{i|\stackrel{}{n}|}{R}})\eta _\stackrel{}{n}+{\displaystyle \underset{l=e,\mu ,\tau }{}}\overline{h}_lL_l\stackrel{~}{\mathrm{\Phi }})+\mathrm{H}.\mathrm{c}.]\}`$
where $`\overline{h}_l=h_l\frac{M_F}{M_\mathrm{P}}`$ , $`|\stackrel{}{n}|^2_{i=1}^\delta n_i^2`$.

As it is shown in Eq. (3), there are small mixings mass terms, $`\overline{h_l}v`$, between any KK state $`\eta _\stackrel{}{n}`$ and left neutrinos. It was firstly noticed in , the four-dimensional Yukawa couplings $`\overline{h}_l`$ are naturally suppressed by the volume factor $`M_F/M_\mathrm{P}`$.

After diagonalizing the mass matrix these terms give rise to the small admixture of the left neutrinos to the heavy KK states, $`\eta _\stackrel{}{n}^{}`$, with the mass $`\sqrt{m^2+\stackrel{}{n}^2/R^2}`$

$$\eta _\stackrel{}{n}^{}\eta _\stackrel{}{n}+B_{l,\stackrel{}{n}}\nu _L^l,B_{l,\stackrel{}{n}}h_l\frac{vM_F}{M_P\sqrt{m^2+\stackrel{}{n}^2/R^2}}.$$
(4)
Applaying to the unitarity condition we arrive to the following light states

$$\nu _{light}^l\frac{1}{\sqrt{1+_\stackrel{}{n}|B_{l,\stackrel{}{n}}|^2}}(\nu _L^l+\underset{\stackrel{}{n}}{}B_{l,\stackrel{}{n}}\eta _\stackrel{}{n}^{}).$$
(5)
The most striking feature of the higher-dimensional scenario, as it follows from Eq.(5), is the loss of the lepton universality in the electroweak processes.

It is useful to define the mixing parameters

$`(s_L^{\nu _l})^2{\displaystyle \underset{\stackrel{}{n}}{}}|B_{l,\stackrel{}{n}}|^2h_l^2{\displaystyle \frac{v^2}{M_F^2}}{\displaystyle \underset{\stackrel{}{n}}{}}{\displaystyle \frac{M_F^4M_P^2}{(m^2+\frac{\stackrel{}{n}^2}{R^2})}}`$
$`\{\begin{array}{c}\frac{\pi h_l^2v^2}{M_F^2}\mathrm{ln}\left({\displaystyle \frac{M_F^2}{m^2}}+1\right),\delta =2\hfill \\ \\ \frac{S_\delta }{\delta 2}\frac{h_l^2v^2}{M_F^2},\delta >2.\hfill \end{array}`$ (9)
In order to evaluate the sumation in the above equation we approximate the sum over the KK states by an integral, which has an upper ultra-violet (UV) cutoff at $`M_FR`$, above which string-threshold effects are expected to become more important.

As can be seen from Table 1, the mixings $`(s_L^{\nu _l})^2`$ may be constrained by a number of new-physics observables induced at the tree level. These observables measure possible non-universality effects in $`\mu `$, $`\tau `$ and $`\pi `$ decays. In this respect, in Table 1 we have defined $`R_\pi =\mathrm{\Gamma }(\pi e\nu )/\mathrm{\Gamma }(\pi \mu \nu )`$ and $`R_{\tau \mu }=B(\tau e\nu \nu )/B(\tau \mu \nu \nu )`$.

Owing to the tower of the KK singlet neutrinos which acts cumulatively in the loops, significant universality-breaking as well as flavour-violating effects are induced in electroweak processes involving $`\gamma `$- and $`Z`$\- boson  interactions. In particular, as has been explicitly shown recently, we find that the cumulative presence of the KK states leads to an effective theory whose Yukawa interactions are mediated by order-unity Yukawa couplings of the original Lagrangian before compactification. In this case, we expect a higher-dimensional non-decoupling phenomenon analogous to the one studied earlier in renormalizable 4-dimensional theories. For example, the effective lepton-flavour-violating vertex $`Zll^{}`$ that occurs in $`\mu eee`$ exhibits the dependence: $`𝒯(Zl^{}l)h_lh_l^{}(v^2/M_F^2)_{k=e,\mu ,\tau }(h_k^2v^2)/M_W^2`$, i.e. its strength increases with the fourth power of the higher-dimensional Yukawa couplings. This should be contrasted with the respective photonic amplitude $`𝒯(\gamma l^{}l)h_lh_l^{}v^2/M_F^2`$, whose strength increases only quadratically.

Based on this cumulative non-decoupling effect, we are able to derive strong limits on the $`M_F/h^2`$, for $`h1`$. As is displayed in Table 2, the strongest limits are obtained from $`\mu \overline{)}eee`$ and the absence of $`\mu `$-to-$`e`$ conversion in nuclei. Of course, the limits presented here contain some degree of uncertainty, which is inherent in all effective non-renormalizable theories with a cut-off scale, such as $`M_F`$. Nevertheless, our results are very useful, since they indicate the generic size of constraints that one has to encounter in model-building considerations with sterile neutrinos.

Now we concentrate in the minimum brane-inspired scheme in which all elements required to explain the neutrino anomalies (the LSND/HDM as well as the solar and atmospheric mass scales) are generated by the physics of extra dimensions . For definiteness, we will be considering a model that minimally extends the standard field content by one bulk neutrino, $`N(x,y)`$ (with zero high dimensional mass term), singlet under the SM gauge group. This propagates on a $`[1+(3+n)]`$-dimensional Minkowski space with $`\delta n`$.

After integrating out heavy Kaluza-Klein states, the effective Lagrangian has the following form

$$_{\mathrm{eff}}=_{\mathrm{SM}}+(\underset{l=e,\mu ,\tau }{}\overline{h}_lL_l\stackrel{~}{\mathrm{\Phi }}\xi _0+\mathrm{H}.\mathrm{c}.),$$
(10)
where $`\xi _0`$ is the zero mode of the Kaluza-Klein states and now

$$\overline{h_l}=\left(\frac{M_F}{M_\mathrm{P}}\right)^{\frac{n}{\delta }}h_l.$$
(11)
In order to account for the LSND or hot dark matter mass scale $`m_\nu 1`$ eV or so, we choose $`\delta =6`$ and $`n=4`$, giving $`\overline{h}_l10^{10}h_l`$.

Now we turn to the Majorana masses for neutrinos. These are crucial in order to generate the mass splittings required in neutrino oscillation interpretations of the solar and atmospheric neutrino anomalies found in underground experiments. As shown in Ref.  the neutrinos may get small Majorana masses via interactions with distant branes where fermion number is maximally broken.

The Majorana part of the neutrino masses is then expected to be

$$m_{ll^{}}f_{ll^{}}\frac{v^2}{M_F}\left(\frac{M_F}{M_\mathrm{P}}\right)^{2n/\delta 4/\delta }$$
(12)
The neutrino mass matrix takes in the basis ($`\nu _e,\nu _\mu ,\nu _\tau ,\nu _s`$ ) ($`\nu _s=\xi _0`$) the form

$$_\nu =\left(\begin{array}{cc}𝓂_{𝓁𝓁^{}}& _𝓁\\ _𝓁^{}^𝒯& \mathcal{0}\end{array}\right).$$
(13)
Here $`M_l=\overline{h}_lv`$. Note that the $`m_s`$ entry in eq. 13 has been omitted since the bulk sector where the sterile neutrino $`\nu _s`$ lives is eight–dimensional.

In the limit that Dirac mass terms ($`M_l`$) are much bigger than Majorana mass terms ($`m_{ll^{}}`$), two of the neutrinos are massless and other two form Dirac state with a mass

$$m_\nu m_{LSND/HDM}=\sqrt{M_e^2+M_\mu ^2+M_\tau ^2}$$
(14)
This state is identified by two angles $`\theta `$ and $`\phi `$ defined as

$$\mathrm{sin}\theta =\frac{M_e}{m_\nu },\mathrm{tan}\phi =\frac{M_\mu }{M_\tau }.$$
(15)
The entries $`m_{ll^{}}`$ only arise due to the breaking lepton number on distant branes. In the case $`\delta =6`$ and $`n=4`$ they are suppressed compared to the Dirac mass terms by the factor $`\frac{v}{M_F}\frac{f_{ll^{}}}{h_l}`$. These terms give masses to the lowest-lying neutrinos and also responsible for splitting Dirac state to two Majorana states. For suitable values of the parameters, these are in the right range to have a solution for solar and atmospheric neutrino deficit. More especifically, from the latest fits one needs

$$\mathrm{\Delta }m_{atm}^23.5\times 10^3eV^2$$
(16)
in order to account for the full set of atmospheric neutrino data.

On the other hand the latest global analysis of solar neutrino data characterized by the best-fit points

$$\mathrm{\Delta }m_{LMA}^23.6\times 10^5eV^2$$
(17)
$$\mathrm{\Delta }m_{SMA}^25\times 10^6eV^2$$
(18)
Now assuming naturalness, namely, that masses and splittings are of the same order $`\mathrm{\Delta }m_{atm}m_{ll^{}}\sqrt{\mathrm{\Delta }m_{}^2}`$ and since $`\mathrm{\Delta }m_{atm}^22m_\nu \mathrm{\Delta }m_{atm}`$, one finds

$$m_\nu 0.8eV\mathrm{for}\mathrm{LMA}$$
(19)
$$m_\nu 0.3eV\mathrm{for}\mathrm{SMA}$$
(20)
characterizing the order of magnitude of the LSND/HDM scale in the LMA/SMA cases.

Of course since clearly the solar mass splitting need not coincide exactly with the lightest state masses, the above estimates are meant to be crude order-of-magnitude estimates only. As a result the LSND/HDM scales both in the LMA and in the SMA case can be larger than estimated above.

If we assume that muonic neutrino coupling to the high dimensional spinor is dominant ($`h_e,h_\tau h_\mu 0.1`$ ) the light sterile neutrino, $`\nu _s`$, combines with $`\nu _\mu `$ and form a Quasi–Dirac state, crucial to account for the hint coming from the LSND experiment, and may also contribute to the hot dark matter of the Universe. Apart from the mass of the Quasi–Dirac state , one has the splittings between its components, as well as the masses of two light active states. The splittings between the heavy states and that characterizing the lighter neutrinos will be associated with the explanations of the atmospheric and solar neutrino anomalies, respectively. The atmospheric neutrino deficit is ascribed to the $`\nu _\mu `$ to $`\nu _s`$ oscilations. The solar neutrino problem could be solved via MSW small or large angle $`\nu _e`$ to $`\nu _\tau `$ solutions. This reproduces exactly the phenomenological features of the model proposed in ref. , providing a complete scenario for the present neutrino anomalies.

Acknowledgements: It is pleasure to thanks the organizers for organizing a very interesting workshop. This work was supported by Spanish DGICYT under grants PB98-0693 and SAB/1998-0136 (A.I.) and by the European Union TMR network ERBFMRXCT960090.


text:
# 𝐽/𝜓 suppression in Pb–Pb collisions at CERN SPS

## 1 Introduction

In ordinary nuclear matter, quarks and gluons are confined inside nucleons. Non-perturbative calculations of Quantum Chromo Dynamics predict that, when temperature exceeds a critical value T$`{}_{c}{}^{}`$ 150-180 MeV, nuclear matter should undergo a phase transition into a state of matter in which quarks and gluons are no more confined into hadrons and behave as free particles. Such a state of matter is named Quark-Gluon Plasma (QGP).

Heavy-ion collisions are a powerful experimental tool to investigate nuclear matter under extreme conditions: the formation of the QGP is expected to occur in these collisions if the critical temperature and energy density required for the phase transition are reached. Several probes have been proposed as signatures of the formation of a deconfined state of matter. In particular Matsui and Satz predicted that the $`J/\psi `$ yield would be suppressed in a deconfined medium due to the Debye screening of the attractive colour force which binds the c and $`\overline{c}`$ quarks together. $`J/\psi `$ suppression is a particularly interesting signature of QGP formation because it probes the state of matter in the earliest stages of the collision, since $`c\overline{c}`$ pairs can only be produced at that time. Moreover the $`J/\psi `$ is a tightly bound state that can not be easily broken by interactions with the hadronic medium and therefore it carries its original message through the different stages of the reacting medium.

At the CERN SPS, the NA38 and NA50 experiments have studied $`J/\psi `$ production using the $`\mu ^+\mu ^{}`$ decay channel with incident proton, oxygen, sulphur and lead ions on several targets.

## 2 Experimental setup and data taking conditions

The $`J/\psi `$ is detected via its $`\mu ^+\mu ^{}`$ decay by a dimuon spectrometer which consists of an air-gap toroidal magnet equipped with two sets of multiwire proportional chambers (4 MWPC upstream the magnet and 4 downstream) for muon tracking purposes and with 6 scintillator hodoscopes to provide the dimuon trigger. The covered rapidity window is $`2.8y_{lab}4.0`$.

The centrality of the collision is estimated event by event using three different detectors. The electromagnetic calorimeter measures the neutral transverse energy ($`E_T`$) produced in the interaction in the pseudorapidity window $`1.1\eta _{lab}2.3`$. The zero degree calorimeter measures the energy $`E_{ZDC}`$, essentially carried by the projectile spectator (non-interacting) nucleons. It covers the pseudorapidity interval $`\eta _{lab}6.3`$. Finally the silicon microstrip detector measures the multiplicity and the angular distribution of charged particles in the acceptance window $`1.5\eta _{lab}3.5`$.

A quartz beam hodoscope (BH), placed about 33 m upstream the target is also used to count the incident lead ions and to reject beam pile-up. The NA50 target is a segmented active Pb target: there can be up to 7 subtargets separated by few centimeters of air. Two quartz blades located off the beam axis on the left and right side of each subtarget identify where the interaction has occurred. For a more detailed description of the detectors see Ref. 4.

## 3 Charmonium production in p–A and light ion interactions

NA38 and NA51 experiments collected an extensive set of measurements of $`J/\psi `$ production using p, O and S beams on several targets. Since charmonium production is a hard process and therefore it scales as A$`\times `$B (where A and B are the projectile and the target mass number respectively), it is possible to define a cross-section per nucleon-nucleon collision as B$`{}_{\mu \mu }{}^{}\sigma _{J/\psi }^{}/(A\times B)`$, where Bμμ is the branching ratio of $`J/\psi `$ into two muons, in order to compare charmonium production in different colliding systems.

The standard analysis method has been described in detail in Ref. 4. The number of J/$`\psi `$ events is determined by fitting the invariant mass spectrum of opposite sign muon pairs above 2.9 GeV/c2. The fit is performed including 5 contributions: $`J/\psi `$, $`\psi ^{}`$, Drell-Yan, open charm (i.e. semileptonic decays of D and $`\overline{D}`$ mesons) and combinatorial background from $`\pi `$ and K decays. An example of invariant mass spectrum is presented in fig. 1

In figure 2 the $`J/\psi `$ cross-section per nucleon-nucleon collision as a function of A$`\times `$B is represented: it can be seen that all the data from p-p up to central S-U collisions show a continuous and monotonic $`J/\psi `$ suppression pattern from the lighter to the heavier interacting nuclei. The systematics of $`J/\psi `$ production from p-p to S-U can be parametrized by the simple law $`\sigma _{AB}=\sigma _0(AB)^\alpha `$ with $`\alpha =0.918\pm 0.015`$. Within the framework of the Glauber model, using a simple first order exponential fit, these data lead to a $`J/\psi `$ absorption cross-section of $`5.9\pm 0.6`$ mb (equivalent to $`6.4\pm 0.8`$ mb for a full Glauber calculation). This result can be understood in terms of ordinary nuclear absorption of a preresonant $`c\overline{c}`$ state meant to become later on, if not destroyed, the fully formed $`J/\psi `$. This normal suppression sets the baseline with respect to which we can compare the pattern of $`J/\psi `$ production in Pb-Pb interactions.

## 4 Anomalous J/$`\psi `$ suppression in Pb–Pb collisions

In fig. 2 it can be seen that the $`J/\psi `$ cross-section measured in Pb–Pb collisions lies $``$ 5 standard deviations below the value expected from the fit to the data from p–p to S–U and this result indicates that there is a new suppression mechanism at work.

In order to perform a study of charmonium production as a function of the collision centrality, it is necessary to find a replacement for A$`\times `$B in the definition of the cross section per nucleon-nucleon collision. Since the Drell-Yan (DY) process provides muons in the same mass range as the $`J/\psi `$ and its cross section is expected to be proportional to the number of elementary nucleon-nucleon collisions without any sizeable nuclear effect, the ratio $`B_{\mu \mu }\sigma _{J/\psi }/\sigma _{DY}`$ can be used to study the centrality dependence of the $`J/\psi `$ suppression with the advantage that in the ratio $`\sigma _{J/\psi }/\sigma _{DY}`$ systematic errors related to detector inefficiencies and flux uncertainties cancel. The transverse energy spectrum is then divided into bins: in each of these bins a fit to the invariant mass spectrum of the muon pairs is performed and the $`J/\psi `$ over DY ratio is calculated.

The disadvantage of this analysis method is that there are statistical fluctuations in the $`\sigma _{J/\psi }/\sigma _{DY}`$ ratio, essentially due to the small statistics of the DY sample. In order to overcome this problem, a new independent analysis method has been developed in which the sample of DY events is replaced by the huge sample of minimum bias (MB) events, collected with a beam trigger which fires every time that a non zero energy deposit is detected by the Zero Degree Calorimeter (even for the most central collisions some energy, from produced particles, falls into the ZDC acceptance).

In this anaysis method the number of $`J/\psi `$ events is obtained without any fitting procedure, simply by counting the number of opposite sign muon pairs in the mass range from 2.9 to 3.3 $`GeV/c^2`$, after combinatorial background subtraction. In order to compare the results of this analysis with the ones obtained with the standard method, the minimum bias reference has to be converted into the Drell-Yan reference. This is done multiplying the measured ET spectrum of the MB events by the ratio of the theoretical shapes of the ET spectra of DY and MB, evaluated in the frame of the Glauber model.

Because of the bigger size of the minimum bias sample, the new analysis is no more affected by statistical fluctuations and a finer $`E_T`$ binning is possible. It is also free from most inefficiencies as it is still computed from a ratio of experimental numbers.

The $`J/\psi `$ over DY ratios obtained with the two independent analysis show a good agreement, as it can be seen in fig. 3 where the complete $`J/\psi `$ suppression pattern obtained with the standard (closed circles) and the minimum bias (open points) analysis for the 1996 and 1998 data samples is plotted as a function of the transverse energy ET.

The solid line visible in the figure represents the “ordinary” nuclear absorption of the $`c\overline{c}`$ pair, which accounts for the $`J/\psi `$ yield measured in lighter colliding systems (from p-p to S-U). It can be observed that for the most peripheral Pb–Pb collisions the $`J/\psi `$ suppression can be accounted for by ordinary nuclear absorption while a clear departure from this trend can be seen at $`E_T`$ 40 GeV, which approximately corresponds to an impact parameter of 8 fm. It is also evident that there is no saturation of $`J/\psi `$ suppression for the most central Pb–Pb collisions: a second drop in the $`J/\psi `$ suppression pattern can be observed at $`E_T`$ 90 GeV.

## 5 Conclusions

The combined results of the NA38 and NA50 experiments clearly indicate a step-wise $`J/\psi `$ suppression pattern with no saturation for the most central collisions.

The experimentally observed suppression is much steeper (when the experimental ET resolution is taken into account) than the results of the predictions of all presently available models of $`J/\psi `$ suppression based on the absorption of the meson by interactions with the surrounding hadronic matter.

On the contrary the observed suppression pattern can be naturally understood in a deconfinement scenario. Since $``$ 30–40% of the $`J/\psi `$’s come from radiative decay of $`\chi _c`$ mesons, the first anomalous step in charmonium suppression can be understood as due to the melting in deconfined matter of the $`\chi _c`$, which is less tightly bound than the $`J/\psi `$. In this scenario, the second drop, observed for more central collisions, would be due to the dissolving of the more tightly bound $`J/\psi `$ meson, which requires energy densities grater than the ones necessary to melt the $`\chi _c`$.

Hence it can be concluded that the $`J/\psi `$ suppression pattern observed in the NA50 data provides significant evidence for deconfinement of quarks and gluons in Pb–Pb collisions.


text:
# A quaternionic generalisation of the Riccati differential equation

## 1 Introduction

The Riccati equation

$$u=pu^2+qu+r,$$

where $`p,q`$ and $`r`$ are functions, has received a great deal of attention since a particular version was first studied by Count Riccati in 1724, owing to both its peculiar properties and the wide range of applications in which it appears. For a survey of the history and classical results on this equation, see for example ,  and . This equation can be reduced to its canonical form ,

$$y+y^2=v,$$
(1)
and this is the form that we will consider.

One of the reasons for which the Riccati equation has so many applications is that it is related to the general second order homogeneous differential equation. In particular, the one-dimensional Schrödinger equation

$$^2uvu=0$$
(2)
where $`v`$ is a function, is related to the (1) by the easily inverted substitution $`y=\frac{u}{u}.`$ This substitution, which as its most spectacular application reduces Burger’s equation to the standard one-dimensional heat equation, is the basis of the well-developed theory of logarithmic derivatives for the integration of nonlinear differential equations . A generalisation of this substitution will be used in this work.

A second relation between the one-dimensional Schrödinger equation and the Riccati equation is as follows. The one-dimensional Schrödinger operator can be factorised in the form

$$^2v(x)=(+y(x))(y(x))$$

if and only if (1) holds.

Among the peculiar properties of the Riccati equation stand out two theorems of Euler, dating from 1760. The first of these  states that if a particular solution $`y_0`$ of the Riccati equatuion is known, the substitution $`y=y_0+z`$ reduces (1) to a Bernoulli equation which in turn is reduced by the substitution $`z=\frac{1}{u}`$ to a first order linear equation. Thus given a particular solution of the Riccati equation, the general solution can be found in two integrations. The second of these theorems  states that given two particular solutions $`y_0,y_1`$ of the Riccati equation, the general solution can be found in the form

$$y=\frac{ky_0\mathrm{exp}(y_0y_1)y_1}{k\mathrm{exp}(y_0y_1)1}$$
(3)
where $`k`$ is a constant. That is, given two particular solutions of (1), the general solution can be found in one integration.

Other interesting properties are those discovered by Picard and Weyr . The first is that given a third particular solution $`y_3`$, the general solution can be found without integrating. That is, an explicit combination of three particular solutions gives the general solution. The second is that given a fourth particular solutions $`y_4`$, the cross ratio

$$\frac{(y_1y_2)(y_3y_4)}{(y_1y_4)(y_3y_2)}$$

is a constant.

This article is concerned with a quaternionic generalisation of the Riccati equation and versions of the above-mentioned theorems of Euler corresponding to this generalisation. Some necessary notation will be introduced in Section LABEL:quats. In Section 3 we propose the quaternionic generalisation of the Riccati equation, which is shown to be a good generalisation for various reasons, including that it is related to the three-dimensional Schrödinger equation

$$\mathrm{}u+vu=0$$
(4)
(where $`\mathrm{}`$ is the three-dimensional laplacian) in the same ways as the Riccati equation is related to (2). In Section 4, we turn our attention to cases in which particular solutions of (8) can be found, some of which differ considerably from the one-dimensional case. Finally, in Section 5, generalisations of Euler’s theorems will be stated and proved.

## 2 Preliminaries

The complex numbers and complex quaterions are denoted by $`,()`$ respectively. The latter consists of elements of the form

$$a=\underset{k=0}{\overset{3}{}}a_ki_k$$

where the $`a_k`$ and the base units $`i_k`$ satisfy the following rules of multiplication

$`i_0^2`$ $`=i_0=i_k^2,i_0i_k=i_0i_k=i_k,k=1,2,3,`$
$`i_1i_2`$ $`=i_2i_1=i_3,i_2i_3=i_3i_2=i_1,i_3i_1=i_1i_3=i_2.`$
The complex unit $`i`$ commutes with the $`i_k.`$ Frequently it is useful to consider a quaternion $`a`$ as being the sum of a scalar and a vector part, denoted respectively

$$a_0:=Sc(a),\stackrel{}{a}:=Vec(a)=\underset{k=1}{\overset{3}{}}a_ki_k.$$

Conjugation is defined as follows,

$$\overline{a}:=a_0\stackrel{}{a}$$

and the modulus is

$$\left|a\right|^2=a\overline{a}=a_0^2+a_1^2+a_2^2+a_3^2$$

so that in particular $`\stackrel{}{a}^2=\left|\stackrel{}{a}\right|^2.`$

Note that in terms of scalars and vectors, the quaternionic product can be written

$$ab=(a_0+\stackrel{}{a})(b_0+\stackrel{}{b})=a_0b_0+a_0\stackrel{}{b}+b_0\stackrel{}{a}\stackrel{}{a},\stackrel{}{b}+\left[\stackrel{}{a}\times \stackrel{}{b}\right]$$

where $`a,b`$ is the standard inner product and $`\left[a\times b\right]`$ the standard cross product in $`^3.`$ In particular

$$\{\stackrel{}{a},\stackrel{}{b}\}=2\stackrel{}{a},\stackrel{}{b}$$
(5)
where $`\{a,b\}=ab+ba`$ is the standard anticommutator.

In what follows functions $`g:\mathbb{\Omega }()`$ will be considered, where $`\mathrm{\Omega }`$ is some domain in $`^3`$. The Moisil-Theodoresco operator $`D`$ is defined on differentiable functions $`g`$ as follows:

$$Dg=\underset{k=1}{\overset{3}{}}i_k_kg$$

where $`_k=\frac{}{x_k}.`$ Due to properties of the quaternionic product, this can be written

$$Dg=div\stackrel{}{g}+gradg_0+rot\stackrel{}{g}.$$

Thus it follows that

$`Sc(D\stackrel{}{g})`$ $`=div\stackrel{}{g},`$
$`Vec(D\stackrel{}{g})`$ $`=rot\stackrel{}{g},`$
and for scalar functions $`u`$

$$Du=gradu.$$

The theorem of Leibnitz for this operator is the following: given a differentiable scalar function $`u`$ and a differentiable quaternionic function $`g,`$

$$D(ug)=D(u)g+uD(g).$$
(6)
The logarithmic derivative (in the sense of Marchenko) of a scalar function $`u`$ such that $`u0`$ in $`\mathrm{\Omega },`$ is defined as

$$\stackrel{ˇ}{}u=u^1Du.$$

The function $`\stackrel{ˇ}{}u`$ is a vector. The derivative is logarithmic in the following sense: given two scalar functions $`u_1,u_2`$ that do not vanish in $`\mathrm{\Omega }`$, formula (6) implies that

$$\stackrel{ˇ}{}(u_1u_2)=\stackrel{ˇ}{}u_1+\stackrel{ˇ}{}u_2.$$

## 3
A three-dimensional generalisation of the Riccati equation

The substitution of the Jackiw-Nohl-Rebbi-’t Hooft ansatz  in the self-duality equation can be written  in the following quaternionic form

$$_tg+Dg+\left|g\right|^2=0$$
(7)
where the subscript $`t`$ denotes differentiation with respect to time. This equation has obvious formal similarities with equation (1). In  the relation of (7) to the Füter operator

$$_t+D$$

was shown. In particular, it was shown that for any $`f\mathrm{ker}(_t+D),`$ the function

$$\frac{2(gradf_0div\stackrel{}{f})}{f_0}$$

is a solution of (7), that is, a class of instantons was obtained. In what follows we will concentrate on the case of time independent, purely vectorial quaternionic functions, but the nonhomogeneous equation will be considered.

The following result generalises to three dimensions the relation, mentioned in the introduction, between the one-dimensional Schrödinger operator and the Riccati differential equation via the logarithmic derivative.

###### Proposition 1

$`\phi `$ is a solution of (4) if and only if $`\stackrel{}{f}:=\stackrel{ˇ}{}\phi `$ is a solution of

$$D\stackrel{}{f}+\stackrel{}{f}^2=v$$
(8)
Proof.

Suppose that there exists a function $`\phi `$ such that $`\stackrel{}{f}=\stackrel{ˇ}{}\phi .`$ Applying (6) gives

$`D\stackrel{}{f}`$ $`={\displaystyle \frac{1}{\phi ^2}}\phi ,\phi {\displaystyle \frac{1}{\phi }}\mathrm{}\phi ,`$
$`\stackrel{}{f}^2`$ $`={\displaystyle \frac{1}{\phi ^2}}\phi ,\phi ,`$
so that $`\frac{1}{\phi }\mathrm{}\phi =v,`$ or equivalently $`\mathrm{}\phi +v\phi =0.`$ Conversely, given a solution $`\phi `$ of (4), $`\stackrel{}{f}=\stackrel{ˇ}{}\phi `$ is a solution of (8).

In  and  it was shown that the three-dimensional Schrödinger operator can be factorised in the following way

$$\mathrm{}vI=(D+M^\stackrel{}{f})(DM^\stackrel{}{f})$$

where $`I`$ is the identity operator and $`M^\stackrel{}{f}q:=q\stackrel{}{f},`$ if and only if equation (8) holds.

Thus the two relations between the Riccati equation (1) and the one-dimensional Schrödinger equation (2) mentioned in the introduction have natural counterparts relating (8) and the three-dimensional Schrödinger equation (4). These relationships suggest that (8) can be considered a good generalisation of (1). It should also be noted that if $`\stackrel{}{f}=f_k(x_k)i_k,`$ then (8) is reduced to

$$_kf_k+f_k^2=v.$$

That is, a one-dimensional solution of (8) is a solution of (1). Equation (8) will be refered to as the Riccati PDE.

Note that the scalar and vector components of equation (8) are respectively

$`div\stackrel{}{f}+\stackrel{}{f}^2`$ $`=v,`$
$`rot\stackrel{}{f}`$ $`=0.`$
The second equation implies that for a simply-connected domain $`\mathrm{\Omega },`$ there exists a scalar function $`\phi `$ such that $`\stackrel{}{f}=grad\phi .`$ Substituting this in the first equation gives

$$\mathrm{}\phi +\phi ,\phi =v.$$
(9)
This equivalence of the Riccati PDE with a scalar elliptic partial differential equation will be used frequently in what follows.

It should also be noted that if the function $`v`$ in (8) is zero, the substitution $`\stackrel{}{f}=\stackrel{ˇ}{}\phi `$ reduces the equation to

$$\mathrm{}\phi =0.$$

Thus the homogeneous Riccati equation can be solved explicitly, and its solutions are of the form $`\stackrel{ˇ}{}\phi `$, where $`\phi \mathrm{ker}\mathrm{}.`$ This generalises the highly restricted class of solutions of the homogeneous equation (1), which are precisely functions of the form

$$\frac{1}{x+c},$$

$`c`$ constant.

## 4
Particular solutions of the equation

In the next section it will be shown that given one solution of the Riccati PDE it can be linearised. Thus in this section we discuss some possibilities for obtaining particular solutions of (8). First we note that if the function $`v`$ is of the form

$$v(x)=v_1(x_1)+v_2(x_2)+v_3(x_3),$$

then assuming that $`\stackrel{}{f}=f_1(x_1)i_1+f_2(x_2)i_2+f_3(x_3)i_3,`$ equation (8) reduces to the system of Riccati ordinary differential equations

$$_kf_k+f_k^2=v_k,k=1,2,3.$$

Thus in this case, a particular solution of (8) can be found if and only if each of the above equations can be solved. Obviously, if any of the $`v_k`$’s are zero, this task is greatly simplified. This situation corresponds to the Schrödinger equation (4) with potential $`v`$ of the form given above, in which case the variables can be separated.

The existence of a large class of solutions of the homogeneous equation, as described in Section 3, motivates the following procedure, which reduces (8) to various scalar differential equations. Substituting the sum of two functions $`\stackrel{}{f}_1,\stackrel{}{f}_2C^1(\mathrm{\Omega })`$ into equation (8) gives the following possible decomposition:

$`D\stackrel{}{f}_1+\stackrel{}{f}_1^2`$ $`=v_1,`$
$`D\stackrel{}{f}_2+\stackrel{}{f}_2^2`$ $`=v_2,`$
$`\{\stackrel{}{f}_1,\stackrel{}{f}_2\}`$ $`=v_3,`$
$`v`$ $`=v_1+v_2+v_{3.}`$
In particular, if $`\stackrel{}{f}_1,\stackrel{}{f}_2`$ are solutions of the homogeneous equation, this is reduced to

$$\{\stackrel{}{f}_1,\stackrel{}{f}_2\}=2\stackrel{ˇ}{}\phi _1,\stackrel{ˇ}{}\phi _2=v,$$
(10)
where $`\phi _1,\phi _2`$ are harmonic functions.

If $`\phi _1=x_1,\stackrel{}{f}_1=\frac{i_1}{x_1},`$ this becomes

$$_1\phi _2=\frac{1}{2}vx_1\phi _2$$

which has solution

$$\phi _2=A(x_2,x_3)\mathrm{exp}(\frac{1}{2}vx_1𝑑x_1)$$

where $`A(x_2,x_3)`$ is an arbitrary function. Thus if $`\phi _2\mathrm{ker}\mathrm{}`$ and $`\phi _2`$ satisfies the above equation, the sum

$$\frac{i_1}{x_1}+\stackrel{ˇ}{}\phi _2$$

is a particular solution of the Riccati PDE.

If instead of choosing $`\phi _1`$ as above $`\phi _1=\phi _2`$ is substituted in (10), the eikonal equation

$$2(\stackrel{ˇ}{}\phi )^2=v$$
(11)
results. Thus for a scalar function $`\phi \mathrm{ker}\mathrm{}_3,`$ which is also a solution of the above eikonal equation, $`\stackrel{}{f}=2\stackrel{ˇ}{}\phi `$ is a solution of (8).

###### Example 2

Let $`\phi `$ be the fundamental solution of the laplacian $`\mathrm{}`$,

$$\phi =\frac{1}{4\pi \left|x\right|}.$$

This function is harmonic and positive in any domain $`\mathrm{\Omega }`$ which does not include the origin. Furthermore it satisfies equation (11) with

$$v=\frac{1}{\left|x\right|^2}.$$

Thus

$$\stackrel{}{f}=\frac{2\stackrel{}{x}}{\left|x\right|^2}$$

is a solution of (8).

## 5

Generalisations of Euler’s theorems on the Riccati equation

We now state and prove the generalisations to the Riccati PDE of Euler’s theorems on the Riccati equation that were mentioned in the introduction. The first of these theorems states that given a particular solution of the Riccati differential equation, the equation can be linearised.

###### Proposition 3

(Generalisation of Euler’s first theorem)

Let $`\stackrel{}{h}=grad\xi `$ be an arbitrary particular solution of (8). Then

$$\stackrel{}{f}=\stackrel{}{g}+\stackrel{}{h}$$
(12)
is also a solution of (8), where $`\stackrel{}{g}=\stackrel{ˇ}{}\mathrm{\Psi }`$ and $`\mathrm{\Psi }`$ is a solution of the equation

$$\mathrm{}\mathrm{\Psi }+2\xi ,\mathrm{\Psi }=0,$$
(13)
or equivalently of

$$div(e^{2\xi }\mathrm{\Psi })=0.$$
(14)
Proof. Substituting (12) in (8) gives

$$D\stackrel{}{g}+\{\stackrel{}{h,}\stackrel{}{g}\}+\stackrel{}{g}^2=0$$

or alternatively, using (5)

$$D\stackrel{}{g}2\stackrel{}{h},\stackrel{}{g}+\stackrel{}{g}^2=0.$$
(15)
Note that, as in (8), the vector part of (15) is $`rot\stackrel{}{g}=0,`$ so that

$$\stackrel{}{g}=grad\mathrm{\Phi }$$

for some function $`\mathrm{\Phi }.`$ If $`\mathrm{\Psi }=e^\mathrm{\Phi },`$ this is equivalent to

$$\stackrel{}{g}=\stackrel{ˇ}{}\mathrm{\Psi }.$$

Equation (15), written in terms of $`\mathrm{\Psi },`$ is

$$\frac{1}{\mathrm{\Psi }^2}(\mathrm{\Psi })^2\frac{1}{\mathrm{\Psi }}\mathrm{}\mathrm{\Psi }\frac{2}{\mathrm{\Psi }}\xi ,\mathrm{\Psi }+\frac{1}{\mathrm{\Psi }^2}(\mathrm{\Psi })^2=0,$$

so that (15) is equivalent to

$$\mathrm{}\mathrm{\Psi }+2\xi ,\mathrm{\Psi }=0.$$

Noting that

$$div(e^{2\xi }\mathrm{\Psi })=2e^{2\xi }\xi ,\mathrm{\Psi }+e^{2\xi }\mathrm{}\mathrm{\Psi }=e^{2\xi }(\mathrm{}\mathrm{\Psi }+2\xi ,\mathrm{\Psi }),$$

this equation can be rewritten in the form

$$div(e^{2\xi }\mathrm{\Psi })=0.$$

Equation (13) is the well-known transport equation, which appears for example in the Ray Method of approximations of solutions to the wave equation, coupled with the eikonal equation. Equation (14) appears in various applications, for example in electrostatics , where $`e^{2\xi }`$ is the dielectric permeability and $`\mathrm{\Psi }`$ is the electric field potential, and as the continuity equation of hydromechanics in the case of a steady flow, where $`e^{2\xi }`$ is the density of the medium .

###### Remark 4

From (14) we have that

$$e^{2\xi }\mathrm{\Psi }=rot\stackrel{}{s}$$

for some vector-valued function $`\stackrel{}{s},`$ or

$$\mathrm{\Psi }=e^{2\xi }rot\stackrel{}{s},$$

where $`\stackrel{}{s}`$ must satisfy the condition

$$rot(e^{2\xi }rot\stackrel{}{s})=0$$

as $`e^{2\xi }rot\stackrel{}{s}`$ must be the gradient of some function.

As mentioned in the introduction,

given two particular solutions $`y_1,y_{2\text{ }}`$of the Riccati ordinary differential equation the general solution can be found in one integration. A natural question is whether this property extends to the Riccati PDE. The form (3) of the general solution found in this case, suggests a similar substitution in (8). This line of reasoning gives the following result.

###### Proposition 5

(Generalisation of Euler’s second theorem)

Let $`\stackrel{}{h}_1=grad\xi _1`$ and $`\stackrel{}{h}_2=grad\xi _2`$ be two particular solutions of (8). Then there exists a scalar function $`\phi `$ such that

$$\stackrel{}{f}=\phi =\frac{\stackrel{}{h}_1w\stackrel{}{h}_2}{w1}$$

is also a solution of (8), where $`w=Ae^{\xi _1\xi _2},A.`$

Proof.

Equation (8) is equivalent to

$$\mathrm{}\phi +\phi ,\phi =v$$
(16)
where $`\stackrel{}{f}=grad\phi .`$ The scalar functions $`\xi _1`$ and $`\xi _2`$ are two solutions of (16). Substituting the expression

$$\frac{\xi _1w\xi _2}{w1}$$
(17)
for $`\phi `$ into this equation, where $`w`$ is a scalar function, gives

$`\mathrm{}\phi `$ $`=div\left({\displaystyle \frac{\xi _1w\xi _2}{w1}}\right)`$
$`={\displaystyle \frac{(w1)(w\mathrm{}\xi _1+w\xi _1\mathrm{}\xi _2)w(w\xi _1\xi _2)}{(w1)^2}}`$
and

$$\phi ,\phi =\frac{(w\xi _1)^22w\xi _1\xi _2+(\xi _2)^2}{(w1)^2}.$$

Simplifying and using the fact that $`\xi _1`$ and $`\xi _{2\text{ }}`$are solutions of (16), the equation is reduced to

$$\mathrm{log}w=\frac{w}{w}=(\xi _1\xi _2),$$

so that

$$w=Ae^{\xi _1\xi _2}$$

for an arbitrary constant $`A.`$ It remains to show that the expression (17) is the gradient of some scalar function $`\phi .`$ This is the case if the rotational of (17) disappears. This is shown as follows, where the identities $`rot\phi =0,\stackrel{}{f}\times \stackrel{}{f}=0`$ are used.

$`rot\left({\displaystyle \frac{\xi _1w\xi _2}{w1}}\right)`$ $`=rot\left({\displaystyle \frac{\xi _1w}{w1}}\right)rot\left({\displaystyle \frac{\xi _2}{w1}}\right)`$
$`=grad{\displaystyle \frac{w}{w1}}\times \xi _1+{\displaystyle \frac{w}{w1}}rot\xi _1`$
$`grad{\displaystyle \frac{1}{w1}}\times \xi _2{\displaystyle \frac{1}{w1}}rot\xi _2`$
$`=grad{\displaystyle \frac{w}{w1}}\times \xi _1grad{\displaystyle \frac{1}{w1}}\times \xi _2`$
$`={\displaystyle \frac{(w1)www}{(w1)^2}}\times \xi _1+{\displaystyle \frac{w}{(w1)^2}}\times \xi _2`$
$`={\displaystyle \frac{w}{(w1)^2}}\times (\xi _2\xi _1)`$
$`={\displaystyle \frac{Ae^{\xi _1\xi _2}}{(w1)^2}}(\xi _2\xi _1)\times (\xi _2\xi _1)`$
$`=0.`$
It must be noted that the new solution gained is not necessarily the general solution, but a larger class of solutions. The following example illustrates this.

###### Example 6

Consider (8) with $`v=1,`$

$$D\stackrel{}{f}+\stackrel{}{f}^2=1.$$

Two solutions of this equation are $`h_1=i_1=gradx_1,h_2=i_2=gradx_2.`$ Applying the above result gives the class of solutions

$$\frac{i_1Ae^{x_1x_2}i_2}{Ae^{x_1x_2}1},A,$$

however a third solution $`h_3=i_3`$ is not included in the above expression as a special case.