text:
# Thermodynamic perturbation theory of the phase behaviour of colloid / interacting polymer mixtures

## I Introduction

The structure, rheology and phase behaviour of sterically stabilized colloidal dispersions are strongly affected by the presence of non adsorbing polymer. Nearly fifty years ago Asakura and Oosawa realized that finite concentrations of polymer coils would induce an effective attraction between colloidal particles, of essentially entropic origin, the so-called depletion interaction AO . Since the initial colloid-polymer Hamiltonian involves only repulsive interactions between all pairs of particles, the polymer-induced effective attraction between colloids, which results from tracing out the polymer degrees of freedom, was referred to as “attraction through repulsion” by A. Vrij. For non-interacting (ideal) polymers, the range of the depletion attraction is independent of polymer concentration, and close to the polymer radius of gyration $`R_g`$, while the depth of the attractive well, when two colloids touch, is proportional to polymer concentration. Consequently, one expects that for sufficiently high concentration, and for not too small size ratio $`q=R_g/R_c`$ (where $`R_c`$ is the radius of the spherical colloids), the effective attraction may drive a depletion-induced phase separation into colloid-rich (“liquid”) and colloid-poor (“gas”) colloidal dispersions, similar to condensation in simple fluids. This phase transition, which is in fact a colloid-polymer demixing transition, was investigated by Gast et al. gast , who first calculated the phase diagram from thermodynamic perturbation theory. Their findings were later confirmed by the Monte-Carlo (MC) simulations of Meijer and Frenkel meijer , and the free volume theory of Lekkerkerker et al. lekk . The predicted phase diagrams agree qualitatively with experimental findings for various colloid / polymer mixtures experiments ; ilett .

More recently the question was raised of how interactions between polymer coils would affect the phase behaviour compared to that of ideal polymers warren ; fuchs ; PRL89 ; aarts ; schmidt . The early theoretical investigations into the problem were made at the two-component level, involving an explicit consideration of the polymer coils. However very recently the depletion-driven effective pair-potential between two colloidal particles in a bath of interacting polymers in good solvent, modelled as self-avoiding walk (SAW) polymers, was calculated by MC simulations JCP117 . A simple analytic form, with coefficients determined by the SAW polymer osmotic equation of state and surface tension, yields excellent agreement with the simulation data over a wide range of polymer to colloid size ratio $`q`$, and polymer concentrations JCP117 .

In this paper we use thermodynamic perturbation theory simplliq  to calculate the phase diagram of mixtures of hard sphere colloids and interacting polymers within the effective one-component representation, whereby the colloidal particles interact via the above-mentioned depletion potential induced by the SAW polymers. The results of these calculations can be directly compared to the predictions of recent MC simulations of the two-component representation of the same system PRL89 , which agree quantitatively with recent experiments ramak . Any discrepancies between the phase diagrams obtained within the effective one-component and two-component representations can then be traced back to more-than-two-body depletion interactions between colloidal particles, which are automatically included in the latter representation, but are of course neglected in the pairwise additive effective one-component picture considered in the present paper.

Thermodynamic perturbation theory was previously applied to mixtures of hard sphere colloids and ideal polymers within the effective one-component representation by Gast et al. gast , and improved by Dijkstra et al. dijkstra . Similar calculations were recently published for mixtures of colloids and star polymers for several functionalities $`f`$ (where $`f`$ is the number of identical arms of the star polymer connected at the centre), again within an effective one-component representation, as well as within the two-component description dzubiella . Star polymers are particularly instructive depletants, since upon varying the functionality, they change continuously from linear polymer ($`f=2`$) to hard-sphere like behaviour ($`f\mathrm{}`$). The phase diagrams obtained in ref. dzubiella  for $`f=2`$ are thus, in principle directly comparable to the results presented in this paper. Such a comparison will be made in section $`4`$, but is only tentative, since the calculations in ref. dzubiella  neglect the polymer concentration-dependence of the effective colloid-polymer and the polymer-polymer interactions which are not negligible PRL89 ; JCP114 ; macromol .

## II One and two-component representations

Consider a system of $`N_c`$ spherical colloidal particles of radius $`R_c`$, in a bath of linear polymers which are in equilibrium with a polymer reservoir of fixed chemical potential $`\mu _p`$. The corresponding semi-grand canonical description is schematically represented in figure 1. The colloidal particles interact via the standard hard sphere potential, while each polymer is made up of $`L`$ monomers or segments; segments from the same or different chains are not allowed to overlap. In good solvent, this excluded volume constraint is the only monomer-monomer interaction and for sufficiently large $`L`$, where chemical details become irrelevant, the interacting polymers may be accurately modelled by self avoiding walks on a three-dimensional lattice. The monomers, moreover, are not allowed to penetrate the hard sphere colloids. At finite colloid and polymer concentrations, such a binary mixture of hard spheres and interacting polymers poses a formidable problem to theoreticians and simulators alike.

One coarse-graining strategy which has proved very successful is to trace out individual monomer degrees of freedom in the “polymers as soft colloids” paradigm, whereby the total interaction between two polymer coils, averaged over all monomer configurations, reduces to an effective (entropic) interaction between their centres of mass, which depends on polymer concentration JCP114 . Similarly, one can trace out the monomer degrees of freedom, to derive a state-dependent effective interaction between the hard sphere colloids and the centres of mass of the polymer coils JCP114 ; macromol . This coarse-graining procedure, which amounts to a reduction of the number of degrees of freedom of each polymer from $`𝒪(L)`$ to $`3`$, leads to a two-component representation of “hard” and effective “soft” colloids, which has been exploited in recent MC simulations to determine the phase diagram of colloid / interacting polymer mixtures for several size ratios $`q`$ PRL89 . However, following the Asakura-Oosawa (AO) strategy for non-interacting polymers, one can carry the coarse-graining procedure one step further, by eliminating the polymer degrees of freedom altogether and determining the resulting depletion interactions between the colloidal particles. If this procedure is carried out in the semi-grand canonical ensemble, the total effective interaction energy between $`N_c`$ colloidal particles for any configuration $`\stackrel{}{r}^{N_c}`$ is:

$$V_{cc}^{eff}(\stackrel{}{r}^{N_c})V_{cc}(\stackrel{}{r}^{N_c})k_BT\mathrm{ln}<\mathrm{exp}[\beta V_{cp}]>_{\mu _p,V,T}(\stackrel{}{r}^{N_c})$$
(1)
where $`V_{cc}`$ is the direct colloid-colloid interaction energy, while $`V_{cp}`$ is the total colloid-polymer interaction; the brackets denote a grand-canonical average over polymer degrees of freedom at fixed polymer chemical potential, volume and temperature, and a given colloid configuration. In the model under consideration, $`V_{cc}`$ and $`V_{cp}`$ are the colloid-colloid and colloid-monomer excluded volume interactions, which are pair-wise additive. The depletion interactions between the $`N_c`$ colloidal particles are given by the second term on the r.h.s. of Eq. (1) and are not, a priori, pair-wise additive. However, for sufficiently low colloid concentration, or small size ratio $`q`$, the pair-wise additive contribution dominates.

### II.1 Depletion potential for interacting polymers

The depletion pair potential between an isolated pair of colloidal spheres has been determined as a function of the centre-to-centre distance, and over a range of interacting polymer concentrations covering the dilute and semi-dilute regimes, in ref. JCP117 . The resulting depletion pair potential is accurately reproduced by the following simple semi-empirical form, inspired by the Derjaguin approximation JCP117 :

$$\beta v_s(r)=\{\begin{array}{cc}\pi R_c\gamma _w(\rho )D_s(\rho )\left(1\frac{x}{D_s(\rho )}\right)^2\hfill & ;0xD_s(\rho ),\hfill \\ 0\mathrm{otherwise}\hfill & \end{array}$$
(2)
$`r`$ is the centre-to-centre separation, while $`x=r2R_c`$ is the surface-to-surface distance between the colloids; $`\gamma _w(\rho )`$ is the polymer surface tension near a planar wall, a function of polymer bulk concentration $`\rho =\rho _b`$ determined in ref. JCP116  for SAW polymers; $`D_s(\rho )`$ is the range of the depletion potential given, according to ref. JCP117 , by:

$$D_s(\rho )=\frac{2\gamma _w(\rho )}{\mathrm{\Pi }(\rho )}\frac{R_{AO}^{eff}(q)}{R_{AO}^{eff}(q0)}$$
(3)
where $`\mathrm{\Pi }(\rho )`$ is the osmotic pressure of the interacting polymers taken from renormalization group calculations oono , namely:

$$\mathrm{\Pi }(\rho )=\rho Z(\alpha \varphi )$$
(4a)
where $`\alpha =2.55`$, $`\varphi =4\pi \rho R_g^3/3`$, $`R_g`$ is the polymer radius of gyration at zero density ($`\rho =0`$), and
$$Z(x)=1+\frac{x}{2}\mathrm{exp}\left[0.309\times [(1\frac{1}{x^2})\mathrm{ln}(1+x)+\frac{1}{x}]\right]$$
(4b)
The second fraction on the r.h.s of (3) accounts for finite curvature of the colloid surface, as estimated approximately from the AO model for ideal polymers JCP116 :

$$\frac{R_{AO}^{eff}(q)}{R_g}=\frac{1}{q}\left(\left(1+\frac{6}{\sqrt{\pi }}q+3q^2\right)^{\frac{1}{3}}1\right),$$
(5)
which reduces to $`2/\sqrt{\pi }`$ in the limit $`q0`$.

### II.2 Depletion potential for ideal polymers

For comparative purposes, we also consider the depletion pair potential induced by ideal polymers. Meijer and Frenkel meijer  showed that to a good approximation this is well represented by the Asakura-Oosawa form AO :

$$\beta v_{id}(r)=\{\begin{array}{cc}\frac{4\pi }{3}\rho \sigma ^3\left\{1\frac{3}{4}\frac{r}{\sigma }+\frac{1}{16}\left(\frac{r}{\sigma }\right)^3\right\}\hfill & ;\mathrm{\hspace{0.33em}2}R_cr\mathrm{\hspace{0.33em}2}\sigma ,\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}$$
(6)
where $`\sigma =R_c+R_{AO}^{eff}`$ is the radius of a sphere around the colloids from which the interpenetrable polymers are excluded, but with an effective $`R_{AO}^{eff}`$ calculated from the insertion free energy of a single colloid in a bath of ideal polymer. The radius is given by eq. (5) JCP116 ; for a hard wall it reduces to $`2/\sqrt{\pi }1.128`$, while it monotonically decreases for increasing size ratio $`q`$ since the polymer can deform around the spherical colloid. For the size ratios considered here the curvature effects are small, on the order of a few $`\%`$ aarts .

The two pair potentials (2) and (6) are always attractive. For interacting polymers the range decreases with increasing polymer concentration $`\rho `$, while the latter is constant for ideal polymers. Furthermore, at any given $`\rho `$, the depth of $`v_s(r)`$ is always less than that of $`v_{id}(r)`$, so that the depletion attraction induced by interacting polymers is weaker than for ideal polymers. Representative examples of the depletion potentials (2) are shown in figure 2.

## III Free energy calculations

Our objective is to draw phase diagrams of colloid / polymer mixtures in the ($`\eta _c`$,$`\varphi _p^r`$) plane, where $`\eta _c=4\pi \rho _cR_c^3/3`$ is the colloid packing fraction, while $`\varphi _p^r=4\pi \rho _p^rR_g^3/3`$ is the ratio of the polymer density in the reservoir ($`r`$) over the overlap density $`\rho _p^{}=3/4\pi R_g^3`$; the latter conventionally separates the dilute ($`\rho _p<\rho _p^{}`$) and the semi-dilute ($`\rho _p>\rho _p^{}`$) regimes. To this end we need to calculate the free energies of the various phases. Within the semi-grand canonical ensemble the required free energy is related to the Helmholtz free energy $`F_c(N_c,N_p,V,T)`$ by a Legendre transformation:

$$F(N_c,\varphi _p^r,V,T)=F(N_c,\mu _p,V,T)=F_c(N_c,N_p,V,T)\mu _pN_p$$
(7)
Since the colloidal particles interact via a hard sphere repulsion and an effective, depletion-driven attraction, the natural way forward is to calculate the free energies of the various phases from thermodynamic perturbation theory, using the well-documented hard sphere fluid as a reference system simplliq . To second order in the high-temperature expansion:

$`F`$ $`=`$ $`F_0+F_1+F_2`$ (8)
$`=`$ $`F_0+W_{N_c}_0+{\displaystyle \frac{\beta }{2}}\left[W_{N_c}^2_0W_{N_c}_0^2\right]`$
where $`F_0=F_0(N_c,V,T)`$ is the free energy of the hard sphere fluid, $`\mathrm{}_0`$ denotes an average over the reference system configurations, and $`W_{N_c}`$ is the perturbation potential energy:

$$W_{N_c}=\underset{i<j}{\overset{N_c}{}}w(r_{ij})$$
(9)
with $`r_{ij}=|\stackrel{}{r}_i\stackrel{}{r}_j|`$, the distance between the centres of colloids $`i`$ and $`j`$; $`w=w(r;\varphi _p^r)`$ is the polymer concentration-dependent depletion potential (2) for interacting polymers ($`w=v_s`$) or (6) for non-interacting polymers ($`w=v_{id}`$). We stress that within the semi-grand canonical description the depletion potential must be calculated for the polymer density in the reservoir, which is unequivocally determined by fixing the chemical potential $`\mu _p`$.

$`F_0`$ in eq. (8) is calculated from the Carnahan and Starling equation of state for the hard sphere fluid CS  while for the hard sphere solid we adopt Hall’s equation of state hall . $`F_1`$ is easily expressed in terms of the hard sphere pair distribution function $`g_0(r)`$ for which we adopted the Verlet-Weis parametrisation in the fluid VW , and the form proposed by Kincaid and Weis for the FCC solid phase KW . The calculation of $`F_2`$ involves three and four-body contributions of the reference system. We have adopted the approximate expression due to Barker and Henderson barkerhend , which only involves the pair distribution function and the compressibility of the reference system. Gathering results:

$$\frac{\beta F}{N_c}=\frac{\beta F_0}{N_c}+\frac{1}{2}\beta \rho _cd^3rg_0(r)w(r)\frac{1}{4}\left(\frac{\rho _c}{p}\right)_0\beta \rho _cd^3rg_0(r)w^2(r).$$
(10)
Note that in the solid phase, $`g_0(r)`$ is the orientationally-averaged pair distribution function.

In order to assess the accuracy of $`F_2`$ and the convergence of the perturbation series (10), we have carried out MC simulations to compute explicitly the fluctuation term in eq. (8), which is approximated by the last term in eq. (10), as well as the total excess free energy. The latter is most conveniently calculated by the standard $`\lambda `$-integration procedure simplliq , whereby the depletion-induced perturbation $`W_{N_c}`$ is gradually switched on, resulting in:

$$F(\lambda =1)=F(\lambda =0)+_0^1W_{N_c}_\lambda 𝑑\lambda $$
(11)
where $`F(\lambda =1)`$ is the required free energy of the fully interacting colloid / polymer mixture, $`F(\lambda =0)F_0`$ is the free energy of the hard sphere reference system, and $`W_{N_c}_\lambda `$ is the statistical average of the perturbation weighted by the Boltzmann factor appropriate for a system of particles interacting via the hard sphere repulsion and the partially switched on depletion potential $`\lambda w(r)`$. The calculation of the free energy $`F`$ hence involves several MC simulations to determine $`W_{N_c}_\lambda `$ for a series of discrete values of $`\lambda `$PR184 , typically $`\lambda =n\times 0.05`$ ($`1n20`$).

The convergence of the perturbation series (10) is illustrated in figure 3 for the depletion potential (2) and a size ratio $`q=0.67`$. As expected, the convergence of the series is faster when the size ratio $`q`$ is larger, the polymer concentration is lower and the colloid packing fraction higher.

The accuracy of each term in the series, and of the truncated sum, were tested by MC simulations of periodic samples of $`N_c=108`$ colloidal particles. Representative results are presented in tables I and II. Table I shows that values of $`F_1`$ from eq. (10) are very close to the simulation data, both for the fluid and the solid. The Barker-Henderson approximation underestimates the absolute value of $`F_2`$, by less than a factor of two in the fluid phase, and by a much larger factor in the FCC crystal phase, where it is totally inadequate. However, as shown in Table II, the sum of the first three terms of the perturbation series (10) yields a total free energy which is surprisingly close to the “exact” MC results from the $`\lambda `$-integration. Similar comparisons show that the predictions of perturbation theory are very reliable both for larger and for somewhat smaller values of $`q`$ (say for $`q0.3`$), but that the predictions rapidly deteriorate for small $`q`$, corresponding to narrow potential wells, as expected. This failure will be illustrated in the case $`q=0.1`$ at the end of the following section.

## IV Phase diagrams

Once the free energies of the fluid and FCC solid phases have been calculated from thermodynamic perturbation theory, as explained in the previous section, the phase diagrams can be calculated using the standard double-tangent construction. Since the initial two-component system is athermal, the depletion potentials (2) and (6) are purely entropic, so that the temperature scales out in the Boltzmann factor, and the resulting phase diagrams are independent of temperature. The phase diagrams for mixtures of colloids and interacting polymers in the ($`\eta _c`$,$`\varphi _p^r`$) plane are shown in figure 4 for four values of the size ratio $`q`$, in the range $`0.1q1`$. The phase diagrams look superficially similar to earlier results obtained for mixtures of colloids and ideal polymers lekk ; dijkstra  or star polymers dzubiella . In particular, for the smaller size ratios, the fluid-fluid phase separation is metastable, and preempted by phase coexistence between a high density solid and a single low density fluid phase. Such a behaviour is a typical signature of “narrow” potential wells like those pictured in figure 2 gast ; lekk . For larger size ratios a stable fluid-fluid phase separation appears with a critical point and a triple point, and the resulting phase diagrams are not unlike those of simple atomic systems in the density-temperature plane (with $`1/\varphi _p^r`$ playing the role of $`T`$).

The corresponding phase diagrams for mixtures of colloids and ideal polymers, calculated using the ideal depletion potential (6) are shown figure 5 for comparison. While they look qualitatively similar to those for interacting polymers in figure 4, there are a number of striking quantitative differences. Because the depletion attraction for ideal polymers (eq. 6) is stronger than that for interacting polymers (eq. (2)) for the same polymer concentration, the fluid-fluid phase separation becomes stable at a larger $`\varphi _p^r`$ for the interacting than for the non-interacting polymers. While the phase diagrams for $`q0.1`$ are fairly close, the differences grow with increasing $`q`$. For $`q1`$ the critical point in figure 5 (ideal case) is at ($`\eta _c=0.18`$,$`\varphi _p^r=0.48`$) compared to (0.25,1.21) in the interacting case, while the triple points are at (0.47,0.92) and (0.43,1.42) respectively, indicating dramatic changes when going from ideal to interacting polymers. Also note that while the critical polymer concentration is practically independent of $`q`$ (for $`q0.35`$) in the ideal case, it shifts to higher values as $`q`$ increases in the interacting case. On the other hand, the critical colloid packing fraction decreases as $`q`$ increases in the non-interacting case, while it is practically constant for interacting polymers.

All these trends are similar to those reported recently in simulations of the two-component description of mixtures of colloids and SAW polymers PRL89 . A detailed comparison between the present perturbation theory results for the effective one-component system, and the phase diagrams determined for the two-component representation is made in figure 6. The agreement between the simulation data for the two-component representation and the predictions of perturbation theory for the effective one-component representation is seen to be reasonable, but not perfect, and to deteriorate as $`q`$ increases. The obvious reason is that perturbation theory only includes the pair-wise additive part of the depletion interactions, while the two-component representation also accounts for effective many-body depletion interactions between colloidal particles. The fact that the phase diagrams obtained from the two-component representation are shifted to lower polymer concentrations relative to the prediction for the effective one-component system indicates that the more-than-two-body depletion interactions are overall attractive in nature<sup>1</sup><sup>1</sup>1Of course part of the difference is also due to the perturbation theory, which, in general, slightly underestimates the value of $`\varphi _p^r`$ along phase boundaries (see e.g. the work of Dijkstra et al.dijkstra ). This suggests that the many-body interactions are slightly more attractive than would be inferred from figure 6. . The opposite trend was found in the case of ideal polymers dijkstra , and is consistent with the more pronounced trends in the limit of large $`q`$ long-pol .

We pointed out earlier that the convergence of thermodynamic perturbation theory is expected to deteriorate when the range of the attractive potential well decreases, i.e. when $`q`$ decreases. To check the reliability of second order perturbation theory at $`q0.1`$, we have systematically computed the “exact” free energy by MC simulations, using the $`\lambda `$-integration (eq. (11)). The phase diagrams determined with the approximate and “exact” free energies are compared in figure 7. The agreement remains acceptable for the fluid-solid transition, even for $`q0.1`$, but the (metastable) critical point of the fluid-fluid transition is at too high a colloid packing fractiondijkstra .

Two-component simulationsPRL89  would be very expensive for small $`q`$, because the number of polymers needed scales as $`q^3`$. However, for sufficiently small $`q`$ we don’t expect many-body interactions to be important, and so our one-component simulation should accurately represent the true colloid / polymer system.

A final instructive comparison is between the present results for the phase behaviour of colloid/interacting polymer mixtures and the results for colloid/star polymer mixtures of functionality $`f=2`$dzubiella , which reduce in fact to interacting linear polymers considered in the present work. The phase diagrams calculated from both depletion potentials within the same approximation (10) for the free energies are compared in figure 8 for similar size ratios $`q`$. Although, as explained earlier in this paper, the depletion pair potential calculated for the $`f=2`$ star polymers dzubiella  is not quite the same as the more accurate one we useJCP117 , the phase-diagrams show the similar trends when compared to ideal polymers.

## V Conclusion

We have shown that traditional thermodynamic perturbation theory, requiring only the well-documented equations of state and pair distribution functions of the fluid and solid phases of the reference hard sphere system, leads to reasonably accurate phase diagrams of mixtures of colloidal particles and interacting polymer coils, provided the appropriate concentration-dependent depletion potential between two colloidal spheres is used. As expected, the agreement between the predictions of the effective one-component description, and the more elaborate two-component description observed for low size ratios $`q`$ deteriorates as $`q`$ increases, due to the enhanced importance of many-body interaction which are neglected in the one-component picture. Nevertheless, the disagreement remains tolerable even at $`q1`$, and in view of the excellent agreement between the predictions of the two-component description PRL89  and recent experimental data ilett ; ramak , we conclude that the effective one-component picture, in conjunction with standard thermodynamic perturbation theory, provides a reliable prediction of the phase diagrams of colloids/polymer mixtures in good solvent.

A direct comparison between the phase diagrams for interacting and ideal polymers calculated at the same level of approximation shows considerable quantitative, and even qualitative differences between the two depletants. The main effect of polymer-polymer interactions is to enhance the miscibility of the colloid / polymer mixtures. Similar conclusions were reached by a number of different recent investigations, based on two-component approaches, including integral equationsfuchs , “polymers as soft colloids”PRL89 , extensions of free-volume theoryaarts , density functional theoryschmidt , and star-polymer potentialsdzubiella . Here we show that the differences between the two types of depletants can be rationalised within a one-component effective potential picture, mainly because for a given $`R_g`$ and density $`\rho _p`$, the depletion potentials for interaction polymers are less attractive than those for interacting polymers.

The results of the present work apply to polymers in good solvent, for which the SAW model constitutes an excellent representation. We plan to examine the situation where solvent quality is such that attractive forces between monomers can no longer be neglected krakovi . Upon lowering the temperature from very high (corresponding to the SAW limit) to the $`\theta `$ temperature, we should be able to investigate the gradual change in the phase diagrams from the fully interacting case to one similar to the ideal polymer limit, which have both been considered in the present paper.

###### Acknowledgements.

The financial support of the EPSRC (grant number RG/R70682/01) is gratefully acknowledged. A.A. Louis acknowledges financial support from the Royal Society. B. Rotenberg acknowledges financial support from the Ecole Normale Supérieure de Paris.


format: text

text:
# A Few Exactly Solvable Models For Riccati Type Equations

## 1 Introduction

The Riccati Equation

$$y^{}=p(t)+q(t)y+r(t)y^2$$
(1)
appears in several branches of applicable mathematics, for example population dynamics and mathematical physics, as in quantum scattering.It is known that it exhibits chaotic behavior .It is also well known that, it reduces to a second order linear differential equation by the substitution, $`y=\frac{u^{}}{u}`$,
when we get,

$$u^{\prime \prime }+g(t)u^{}+f(t)u=0$$
(2)
(1) and (2) have been studied in detail for a long time.
It may be mentioned that if two solutions of the Riccati equation are equal at a point, then they coincide.The reason lies in the continuity of the logarithmic derivative $`\frac{u^{}}{u}`$ given in the substitution that lead from (1) to (2).
We will now consider exact solutions of the Riccati equation in two different contexts,from the field of population dynamics and the field of quantum scattering.

## 2 Some exact solutions

### 2.1 A Problem In Population Dynamics

We now examine coupled Riccati type equations of relevance to a population model.

#### Gause’s Model

Consider two-species populations occurring together,and assume that the growth of each is inhibited by members ,both of its own and of the other species.Denoting the number of individuals in species 1 as $`N_1`$ and species 2 as $`N_2`$,we have the Gause’s competition equations:

$$\frac{1}{N_1}\frac{dN_1}{dt}=r_1a_1N_1a_2N_2$$
(3)
$$\frac{1}{N_2}\frac{dN_2}{dt}=r_2a_3N_1a_4N_2$$
(4)
where $`r_1`$,$`r_2`$,$`a_1`$,$`a_2`$ are defined below. Thus we are assuming that the per capita growth of each population at an instant is a linear function of the sizes of the two competing populations at that instant.Each population would grow logistically if it were alone with logistic parameters $`r_1`$ and $`a_1`$ for species 1 and $`r_2`$ and $`a_2`$ for species 2 .
In general the simultaneous differential equations cannot be explicitly solved.We now consider a particular set of circumstances in which they can be solved.Thus for (3) and (4) we specialize to, with a more convenient notation ,

$$\frac{dy}{dt}=1xy^2$$
(5)
$$\frac{dx}{dt}=1+yx^2$$
(6)
From (5) and (6) we get,

$$\frac{dy}{dx}=\frac{1xy^2}{1+yx^2}$$
(7)
Integrating (7) we have

$$yx=x^2y^2+C$$
(8)
From (8) when $`y=0`$,$`C=x_0`$ so that $`y<x`$
Reverting back to $`N_1`$ and $`N_2`$, this is

$$N_1N_2=N_{2}^{}{}_{}{}^{2}N_{1}^{}{}_{}{}^{2}+C$$
(9)
From (9) when $`N_1=0`$,$`C=N_0`$ so that $`N_1<N_2`$.
The Figure below illustrates a particular case of the above solution with the two populations $`N_1`$ and $`N_2`$ along Y and X axes.

### 2.2 Problems from Scattering Theory

By reversing the method used above ,the quantum mechanical radial Schrödinger equation

$$u^{\prime \prime }H(r)u=0$$
(10)
can be reduced to a Riccati equation

$$v^{}+v^2=H(r)$$
(11)
by the substitution $`v=\frac{u^{}}{u}`$. The form (10)or (11) is used in phase shift analysis for example in Calogero’s variable phase approach .It is also possible to use the form (11) for building up an iterative procedure .We would now like to point out that (11) can be used in a different context for providing exact solutions for specified classes of the potential function f(r). Let us write in (11)

$$v\frac{1}{f}$$

$$f^{}=1gf$$

Whence we get

$$H=\frac{g}{f}=\frac{ge^{{\scriptscriptstyle g𝑑r}}}{e^{{\scriptscriptstyle g𝑑r}}+C}$$
(12)
1.We put g=$`\alpha `$ in (12)

$$H=\frac{\alpha ^2}{1+c\alpha e^{\alpha r}}$$

The potential H is now of the form

$$H=\frac{D}{A+Be^{\alpha r}}$$

This is the well-known Wood-Saxon potential. The solution u is given by

$$u=C_1(e^{\alpha r}+C_2\alpha )$$

2.Putting g=r in (12), the potential H is given by

$$H=\frac{re^{r^2/2}}{e^{r^2/2}𝑑r+c_1}$$

This is a modified Gaussian potential and the solution u is given by

$$u=Ke^{r^2}𝑑r+C_3$$

The Graph is shown in the figure below and is seen to fall very steeply indicative of a confined state or particle.

3.Putting g=1/r in (11),the potential H is given by

$$H=\frac{2}{r^2+c_4}$$

This is a shifted inverse square potential and the solution u is given by

$$u=K^{}(r^2+C_4)$$

## 3 Other Exact and Asymptotic Solutions

We now consider some asymptotic solutions of (11),which we write as,

$$u^{}=u^2+f$$
(13)
Let

$$f𝑑rg,$$

f and g being bounded functions. So (13) can be written as ,

$$z^{}=(z+g)^2>0,z=ug,$$
(14)
which shows that z is an increasing function of r.Suppose z is unbounded.So for large r,we should have from (14)

$$zz^2$$

whence

$$z=\frac{1}{cr}$$

which $`0`$ as $`r\mathrm{}`$.
This is a contradiction. Therefore z is bounded and so also u ,that is $`uM`$ as $`r\mathrm{}`$
Therefore $`u^{}0`$ (asymptotically)
Therefore for large r,(13) becomes

$$u^2=f,\text{whence}$$

$$u=\pm (f)^{1/2}$$
(15)
By way of verification of (15),let us consider (13) with,

$$f=\frac{\alpha ^2}{r^2}$$

So we expect that for large r, $`u`$ $`\pm `$ $`\frac{\alpha }{r}`$ by (15). Let us put

$$u=\frac{\beta }{r}$$
(16)
in (13). So

$$\frac{\beta }{r^2}=\frac{\beta ^2}{r^2}\frac{\alpha 2}{r^2}$$

or

$$\beta =\frac{1\pm \sqrt{1+4\alpha ^2}}{2}\pm \alpha \text{if}\alpha 1$$

That is

$$u=\pm \frac{\alpha }{r}$$

(everywhere, and so also for large r). Using this example,with transformations of the independent variable,we can generate similar solutions. For example, if we substitute for r, $`t=t(r)`$, (13) becomes

$$t^{}\dot{u}=u^2+f,$$

where the dot denotes the derivative with respect to t.The choice $`t^{}=f(r)`$, leads to a similar equation, and one can verify that for

$$f(r)=\frac{a}{3}r^{2/3}a^2r^{2/3},u=ar^{1/3}$$

is a solution.More generally as can be easily verified $`ar^n`$ is a solution for $`f(r)=nar^{n1}a^2r^{2n}`$ and so on.However in these examples, neither u nor u’ are asymptotically bounded.
Finally it maybe observed that if one solution of the Riccati equation (13) is known,then others could be derived therefrom.


text:
# Neutron-gravity interferometry experiment: testing Earth’s rotating frequency fluctuations via neutron Berry’s phases due to spin-rotation couplings

## Abstract

$`𝒜BSTRACT`$

Neutron spin can be coupled to the Earth’s rotating frequency. Once if the Earth’s rotating frequency is time-dependent, then the neutron will acquire a Berry’s topological phase (cyclic adiabatic geometric phase). So, a potential method to investigate the Earth’s time-varying rotating frequency by measuring the phase difference between geometric phases of neutron spin polarized vertically down and up is proposed.

In gravity theory, the intrinsic spin of a particle can be regarded as a gravitomagnetic dipole moment. For this reason, it can be coupled to the gravitomagnetic fields. According to the principle of equivalence, the rotating frequency of a non-inertial system can be considered a gravitomagnetic field, and in consequence the Coriolis force acting upon a moving particle can be regarded as a Lorentz (gravitomagnetic) force. Thus, a spinning particle in a rotating frame of reference will be coupled to the rotating frequency, which is termed the spin-rotation coupling. The Hamiltonian describing this interaction of particle spin $`𝐒`$ with the rotating frequency $`\stackrel{}{\omega }`$ of a non-inertial frame is written as $`H=\stackrel{}{\omega }𝐒`$. The spin-rotation coupling leads to the inertial effects of the intrinsic spin of a particle, for instance, although the equivalence principle still holds, the universality of Galileo’s law of freely falling particles is violated, provided that the spin of falling particle is polarized vertically up or down in the non-inertial frame.

If the Hamiltonian $`H`$ of a quantum system is time-dependent, then it will give rise to a geometric phase (Berry’s phase). Differing from dynamical phase that depends on dynamical quantities of systems such as energy, frequency, velocity as well as coupling coefficients, geometric phase is independent of these dynamical quantities. Instead, it is only related to the geometric nature of the pathway along which quantum systems evolve. This, therefore, implies that geometric phase presents the topological and global properties of quantum systems in time-evolution process, and that it possesses the physical significance and can thus be applied to various fields of physics.

If the Earth’s rotating frequency is time-dependent, i.e., $`\stackrel{}{\omega }(t)`$ whose precession frequency and precessional cone angle are respectively $`\mathrm{\Omega }`$ and $`\theta `$, then the interacting Hamiltonian of neutron spin-rotation coupling is $`H(t)=\stackrel{}{\omega }(t)𝐒`$. The time-dependent quantum system is governed by the time-dependent Schrödinger equation. By using Berry’s phase formula or Lewis-Riesenfeld invariant theory, one can arrive at the expression for Berry’s phase, $`\varphi _\pm ^{(\mathrm{Berry})}(T)=\pm \frac{1}{2}\mathrm{\Omega }(1\mathrm{cos}\theta )`$ in one precessional cycle ($`T=\frac{2\pi }{\mathrm{\Omega }}`$), acquired by the neutrons, where the sign $`\pm `$ corresponds to the cases of neutron polarized vertically up and down, respectively.

The variation of the Earth’s rotating frequency may be caused by the motion of interior matter, tidal force and the motion of atmosphere as well. Once we have information on the Earth’s rotating frequency, it is possible to investigate the motion of matter on the Earth. Here we suggest a potential approach to detecting the fluctuations of the Earth’s time-dependent rotation: specifically, measuring the geometric phase arising from the interaction between neutron spin and the Earth’s rotation by using neutron interferometry experiment, which has been successful in investigating the Aharonov-Carmi effect, the gravitational analogue to the Aharonov-Bohm effect in electrodynamics.


text:
# Transmission resonances and supercritical states in a one dimensional cusp potential

## I Introduction

The study of scattering and bound states of nonrelativistic particles in the Schrödinger equation framework is a well known and understood problem Newton . The scattering of low momentum particles by one-dimensional potentials that are well behaved at infinity has zero transmission coefficient at zero energy, and the reflection coefficient is unity unless the potential supports a half bound state (zero energy resonance). In that case there is no reflection and the transmission coefficient is unity provided that the potential is symmetric. The pioneering works on the computation of transmission and reflection coefficients in one-dimensional potentials go back to Bohm Bohm  who coined the definition of transmission resonance to scattering states with transmission coefficient equal to unity and vanishing reflection. Recently Senn ; Sassoli , these results have been generalized to asymmetric potentials.

The generalization of the concept of transition resonances to relativistic quantum mechanics presents some subtleties related to the Dirac equation Kennedy2 ; Dombey . For a free Dirac particle there are positive energy states $`(E>m)`$ as well as negative energy states $`(E<m)`$. In the presence of an external potential, half bound states occur for $`E=m`$ and $`E=m`$, in contrast to the non-relativistic case where these only exist at $`E=0.`$ In the Dirac equation framework one should talk about zero momentum resonances rather than zero energy resonances.

The study of scattering and bound states in the presence of strong electromagnetic fields in relativistic quantum mechanics is also a rather old problem. Soon after the formulation of the Dirac equation Dirac , Klein Klein  found that for very high potential barriers an unusually large number of electrons penetrate into the wall. Moreover, those electrons have negative energy. Sauter Sauter  found the same results in the more general case where the barrier smoothly grows. Pair creation for the Klein paradox is stimulated by the incoming electron beam. When the depth of a square well exceeds $`2mc^2`$, the energy of the lowest bound state becomes equal to $`mc^2`$ , and the bound level joins the states of the Dirac sea. The bound states dive into the energy continuum Greiner .

The Klein paradox in a one-dimensional Saxon-Woods potential was studied by Dosch et al Dosch , showing that in this case transmission resonances exist. Kennedy Kennedy  discussed supercriticality in this potential, demonstrating that, as in the square potential case, the one-dimensional Woods-Saxon potential has transmission resonances of zero momentum when it supports a half bound state at $`E=m`$.

The Saxon-Woods potential presents the particularity of being a smoothed form of a square well; the Dirac equation presents half bound states with the same asymptotic behavior as those obtained with a potential barrier Kennedy .

In the present article we discuss the problem of computing transmission resonances and supercriticality of the one-dimensional symmetric screened potential (Fig 1):

$$U(x)=V_0\mathrm{exp}(\left|x\right|/a).$$
(1)
The potential (1) is an asymptotically vanishing potential for large values of the space variable $`x`$. The parameter $`V_0`$ shows the height of the barrier or the depth of the well depending on its sign. The positive constant $`a`$ determines the shape of the potential. A noticeable difference between the Woods-Saxon potential and (1) is that this one does not exhibit a square barrier limit. The potential (1), for $`V_0<0,`$ can be regarded as a screened one dimensional Coulomb potential Dominguezadame .

It is the purpose of the present article to show that the potential (1), for $`V_0>0`$, supports transmission resonances. We solve the 1+1 Dirac equation explicitly and compute the bound energy states for $`U(x),`$ for $`V_0<0,`$ and show that this potential supports supercritical states which correspond to transmission resonances for the potential (1), for $`V_0>0.`$

The article is structured as follows. In Sec. II we solve the 1+1 Dirac equation in the potential (1). In Sec III we compute the bound states and the condition that the potential should satisfy in order to support supercritical states $`(E=m)`$. In Sec. IV. we compute the scattering states and show that for $`k=0`$ the condition for transmission resonances corresponds to that of supercriticality for $`V_0<0`$.

## II The Dirac equation

The Dirac equation in the presence of an external electromagnetic field $`A_\mu `$ has the form

$$\left\{\gamma ^\mu (_\mu ieA_\mu )+m\right\}\mathrm{\Psi }=0,$$
(2)
where the four-vector potential (1) can be written in a covariant way as

$$A^\mu =V_0\mathrm{exp}(\left|x\right|/a)\delta _0^\mu .$$
(3)
In two dimensions, the Dirac equation (2) in the presence of the spatially dependent electric field (3) reads

$$\left\{\gamma ^0(_t+ieV_0\mathrm{exp}(\left|x\right|/a))+\gamma ^1_x+m\right\}\mathrm{\Psi }=0$$
(4)
where the $`\gamma ^\mu `$ matrices satisfy the commutation relation $`\{\gamma ^\mu ,\gamma ^\nu \}=2\eta ^{\mu \nu }`$ and the metric $`\eta ^{\mu \nu }`$ has the signature $`(,+).`$ Here and throughout the present article we adopt units where $`\mathrm{}=c=1.`$

Since we are working in 1+1 dimensions, it is possible to choose the following representation of the Dirac matrices

$$\gamma ^0=i\sigma ^2=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\gamma ^1=\sigma ^1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$$
(5)
Taking into account that the vector potential (3) does not depend explicitly on time, we can separate the time dependence of the spinor $`\mathrm{\Psi }`$ as follows:

$$\mathrm{\Psi }=\left(\begin{array}{c}\mathrm{\Psi }_1(x)\\ \mathrm{\Psi }_2(x)\end{array}\right)\mathrm{exp}(iEt)$$
(6)
Substituting the gamma matrices (5) into the Dirac equation, (4) we obtain the following system of coupled equations

$$\left(\frac{d}{dx}i(eV_0\mathrm{exp}(\left|x\right|/a)E)\right)\mathrm{\Psi }_1+m\mathrm{\Psi }_2=0$$
(7)
$$\left(\frac{d}{dx}+i(eV_0\mathrm{exp}(\left|x\right|/a)E)\right)\mathrm{\Psi }_2+m\mathrm{\Psi }_1=0$$
(8)
Introducing the new variable

$$y=2iaeV_0\mathrm{exp}(x/a)x>0$$
(9)
we reduce the system (7)-(8) to

$$\left(y\frac{d}{dy}\frac{y}{2}+iEa\right)\mathrm{\Psi }_2=ma\mathrm{\Psi }_1$$
(10)
$$\left(y\frac{d}{dy}+\frac{y}{2}iEa\right)\mathrm{\Psi }_1=ma\mathrm{\Psi }_2$$
(11)
The regular solutions of the system of equations (10)-(11),can be expressed in terms of the Whittaker function $`M_{k,\mu }(y)`$ Abramowitz  as follows

$$\mathrm{\Psi }_2=b_1y^{1/2}M_{k,\mu }(y)$$
(12)
$$\mathrm{\Psi }_1=\frac{1/2+\mu +k}{ma}b_1y^{1/2}M_{k+1,\mu }(y)$$
(13)
where $`b_1`$ is a normalization constant, and

$$k=iEa1/2,\mu =i\sqrt{E^2m^2}a$$
(14)
For negative values of $`x`$, we introduce the auxiliary variable $`\overline{y}`$

$$\overline{y}=2aieV_0\mathrm{exp}(x/a)x<0.$$
(15)
In terms of the new variable $`\overline{y},`$ the system of equations (7)-(8) takes the form

$$\left(\overline{y}\frac{d}{d\overline{y}}+\frac{\overline{y}}{2}iEa\right)\mathrm{\Psi }_2+ma\mathrm{\Psi }_1=0$$
(16)
$$\left(\overline{y}\frac{d}{d\overline{y}}\frac{\overline{y}}{2}+iEa\right)\mathrm{\Psi }_1+ma\mathrm{\Psi }_2=0$$
(17)
and its solution can be written as

$$\mathrm{\Psi }_1=c_1\overline{y}^{1/2}M_{k,\mu }(\overline{y})$$
(18)
$$\mathrm{\Psi }_2=\frac{1/2+\mu +k}{ma}c_1\overline{y}^{1/2}M_{k+1,\mu }(\overline{y})$$
(19)
where $`c_1`$ is a normalization constant

## III Bound states

The potential (3) supports bound states provided that $`m^2E^2>0.`$ In this case one has that $`\mu `$ can be written as

$$\mu =a\sqrt{m^2E^2}$$
(20)
and the condition for existence of bound energy states can be obtained after imposing the normalizability of the wave function as well as the continuity of $`\mathrm{\Psi }`$ at $`x=0.`$ Noticing the normalizability of the spinor $`\mathrm{\Psi }`$ upon the scalar product

$$\mathrm{\Psi },\mathrm{\Phi }=\overline{\mathrm{\Psi }}\gamma ^0\mathrm{\Phi }𝑑x$$
(21)
imposes that the spinor components $`\mathrm{\Psi }_1`$ and $`\mathrm{\Psi }_2`$ should vanish as $`\left|x\right|\mathrm{}`$. This condition is equivalent to imposing the regularity of the spinor at $`y=0,`$ and $`\overline{y}=0.`$ The solutions (12)-(13), and (18)-(19) satisfy this boundary condition. The energy levels can be computed by imposing the continuity of the spinor $`\mathrm{\Psi }`$ at $`x=0.`$

$$c_1M_{k,\mu }(2ieaV_0)=b_1\frac{1/2+\mu +k}{ma}M_{k+1,\mu }(2ieaV_0)$$
(22)
$$b_1M_{k,\mu }(2ieaV_0)=c_1\frac{1/2+\mu +k}{ma}M_{k+1,\mu }(2ieaV_0)$$
(23)
The existence of non-trivial $`c_1`$ and $`b_1`$ satisfying the above system (22)-(23) is equivalent to imposing the condition

$$\left(\frac{1}{2}+\mu +k\right)^2M_{k+1,\mu }^2(2ieaV_0)+m^2a^2M_{k,\mu }^2(2ieaV_0)=0.$$
(24)
The expression (24) shows that the screened potential (1) is able to bind particles, a property that is absent in the purely one-dimensional Coulomb potential. The result (24) is equivalent to that obtained in Ref. Dominguezadame .

As the potential depth increases for large values of $`V_0`$, the energy eigenvalue of any bound state also decreases. When this energy reaches $`E=m`$, the bound states merge with the negative energy continuum and the potential becomes supercritical. When $`E=m,`$ $`\mu =0`$ and we obtain from Eq. (24) the condition for supercriticality:

$$M_{k,0}(2ieaV_0)M_{k,0}(2ieaV_0)M_{k+1,0}(2ieaV_0)M_{k+1,0}(2ieaV_0)=0$$
(25)
The condition (25) permits one to compute the value of $`V_0`$ for a given shape parameter $`a`$. Fig. 2. shows how depends the potential depth on the shape parameter $`a.`$

## IV Scattering states

In order to discuss the transmission and reflection across the barrier given by the potential (3) , for $`V_0>0,`$we identify the incoming traveling wave. Since the vector potential (1) vanishes as $`x\mathrm{}`$, we have the result that the incoming wave should resemble a plane wave for large negative values of $`x`$. Looking at the asymptotic behavior of the Whittaker function $`M_{k,\mu }(z)`$ as $`z0`$

$$M_{k,\mu }(z)\mathrm{exp}(z/2)z^{\frac{1}{2}+\mu },$$
(26)
it is not difficult to obtain that the incoming wave reads

$$\mathrm{\Psi }_{inc}(\overline{y})=\left(\begin{array}{c}\overline{y}^{1/2}M_{k,\mu }(\overline{y})\\ \frac{\frac{1}{2}+\mu +k}{ma}\overline{y}^{1/2}M_{k+1,\mu }(\overline{y})\end{array}\right)$$
(27)
With the help of the relation (26), we obtain that, as $`x\mathrm{},`$ the spinor $`\mathrm{\Psi }_{inc}`$ takes the form

$$\mathrm{\Psi }_{inc}\left(\begin{array}{c}1\\ \frac{\frac{1}{2}+\mu +k}{ma}\end{array}\right)(2ieaV_0)^\mu \mathrm{exp}(ix\sqrt{E^2m^2})$$
(28)
and consequently the spinor (27) satisfies the required asymptotic behavior. Analogously, we have that the reflected wave can be written as

$$\mathrm{\Psi }_{refl}(\overline{y})=\left(\begin{array}{c}\overline{y}^{1/2}M_{k,\mu }(\overline{y})\\ \frac{\frac{1}{2}\mu +k}{ma}\overline{y}^{1/2}M_{k+1,\mu }(\overline{y})\end{array}\right)$$
(29)
as $`x\mathrm{},`$ we have that the reflected wave $`\mathrm{\Psi }_{refl}`$ has the asymptotic form

$$\mathrm{\Psi }_{refl}\left(\begin{array}{c}1\\ \frac{\frac{1}{2}\mu +k}{ma}\end{array}\right)(2ieaV_0)^\mu \mathrm{exp}(ix\sqrt{E^2m^2})$$
(30)
The transmitted wave $`\mathrm{\Psi }_{trans}`$ has the form

$$\mathrm{\Psi }_{trans}(y)=\left(\begin{array}{c}\frac{\frac{1}{2}\mu +k}{ma}y^{1/2}M_{k+1,\mu }(y)\\ y^{1/2}M_{k,\mu }(y)\end{array}\right)$$
(31)
whose asymptotic behavior for large values of $`x`$ is

$$\mathrm{\Psi }_{trans}\left(\begin{array}{c}\frac{\frac{1}{2}\mu +k}{ma}\\ 1\end{array}\right)(2ieaV_0)^\mu \mathrm{exp}(ix\sqrt{E^2m^2})$$
(32)
The computation of the reflection and transmission coefficients can be achieved with the help of the continuity of the spinor wave-function at the boundary $`x=0`$.

$$A\mathrm{\Psi }_{inc}(x=0)+B\mathrm{\Psi }_{refl}(x=0)=C\mathrm{\Psi }_{trans}(x=0)$$
(33)
Substituting (27), (29) and (31) into (33), we obtain

$$\frac{B}{A}=\frac{M_{k,\mu }(2ieaV_0)M_{k,\mu }(2ieaV_0)M_{k+1,\mu }(2ieaV_0)M_{k+1,\mu }(2ieaV_0)}{\frac{(\frac{1}{2}\mu +k)^2}{m^2a^2}M_{k+1,\mu }^2(2ieaV_0)M_{k,\mu }^2(2ieaV_0)}$$
(34)
The condition $`B/A=0`$ is equivalent to considering that there is no reflected wave, i.e, a transmission resonance condition. From Eq. (34) we obtain that a half bound state or zero momentum resonance condition $`\mu =0`$ satisfies the equation

$$M_{k,0}(2ieaV_0)M_{k,0}(2ieaV_0)M_{k+1,0}(2ieaV_0)M_{k+1,0}(2ieaV_0)=0$$
(35)
Equation (35) gives the values of $`V_0`$ and $`a`$ for which zero momentum resonances exist. This condition, with $`E=m,`$ and $`V_0>0,`$ is equivalent to (25) for $`E=m`$ and $`V<0.`$ In fact, after taking the complex conjugate of Eq. (35), using the the relation Gradshtein

$$z^{\frac{1}{2}}M_{k,0}(z)=(z)^{\frac{1}{2}}M_{k,0}(z)$$
(36)
for the Whittaker functions, and substituting the value of $`k`$ defined by the relation (14), we reobtain the supercriticality condition (25). The dependence of $`V_0`$ on $`a`$ for transmission resonance states is shown in Fig. 2.

## V Concluding remarks

The relations (25) and (35) show that the one dimensional screened potential (1) supports supercritical states and consequently half bound states. This result is not trivial in view of the fact that the one-dimensional Coulomb potential presents only scattering states Dominguez2 . The potential (1) does not exhibit a square barrier limit, and for very small values of $`a`$ and constant value of $`V_0a`$ it can be regarded as a delta potential Dominguezadame ; therefore the relations (25) and (35) also hold for a delta potential. The results obtained in this paper show that very cusp symmetric potentials also support half bound states.

.

###### Acknowledgements.

One of the author (VMV) wants to express his gratitude to the Alexander von Humboldt Foundation for financial support.


text:
# Structure of Extremely Nanosized and Confined In-O Species in Ordered Porous Materials

## Abstract

Perturbed-angular correlation, x-ray absorption, and small-angle x-ray scattering spectroscopies were suitably combined to elucidate the local structure of highly diluted and dispersed InOx species confined in porous of ZSM5 zeolite. This novel approach allow us to determined the structure of extremely nanosized In-O species exchanged inside the 10-atom-ring channel of the zeolite, and to quantify the amount of In2O3 crystallites deposited onto the external zeolite surface.

The study of structural, magnetic and electronic properties of nanostructured, subnanostructured and, in the other extreme, highly diluted monoionic species in solids deserves increasing attention not only from a fundamental point of view but also for their technological applications Li . In one extreme, highly-dispersed metallic exchanged-atoms in an infinite variety of compounds present, e.g., very important catalytic properties Greegor , which are not clearly correlated with the active species responsible for them since their physical entities are often unknown. The design of new catalysts with improved activity, selectivity and stability requires the complete knowledge of the local environment of the active centers and its correlation with the desired reaction. The complete characterization of this kind of structures is at present a challenging problem not only in catalysis but also in many fields of physics.

An additional difficulty is the confinement of the diluted species inside porous materials. In general, experimental techniques based on energetic probes that strongly interact with the materials, may destroy them. On the other side, low energetic probes are unuseful since they have to pass through the walls around the “hidden” locations of the confined clusters or atoms and their kinetic energy is attenuated.

The importance of the extended X-ray absorption fine structure (EXAFS) spectroscopy to study structural properties in crystalline solids Greegor , nanoclusters Bazin  and highly dispersed ionic species in catalysis Via  has been long acknowledged. EXAFS constitutes a powerful “atom selective” technique to extract direct information about type, number and distances of neighbors of the absorber atom Teo ; Groot . However, in EXAFS analysis, if more than one species has the same type of bonds (same element and similar bond-lengths), the information from each one would be almost impossible to extract directly, unless some additional information is known.

The perturbed-angular-correlation (PAC) technique requires a suitable probe-atom (native or foreign) to be in the system under study, being in this sense also ”atom selective”. PAC enables the precise determination - at the probe site - of the electric-field gradient tensor (EFG), which is extremely sensitive to the anisotropy of the electronic density near the nucleus, which in turn reflects the probe-neighboring coordination. PAC has been intensively applied to many fields in science during the past two decades with success Lerf  and, very recently, to the development of accurate ab initio calculations of the EFG at diluted impurity sites in crystalline solids Lany . However, structural information is not easily obtained if we are dealing with highly dispersed, disordered, and coexistent multiple-species.

In this letter, we will show that the combination of gamma- and X-ray-based techniques like PAC, EXAFS, and small-angle X-ray scattering (SAXS) can be a novel and powerful experimental tool that can be applied to the emerging field of structure determination of extremely nanosized and confined species, like In-species at exchange sites in ZSM5 zeolites.

The crystalline zeolite framework, an aluminum-silicate structure crossed by channels and cavities, is a challenging laboratory to study the exchange and deposition of confined particles in internal surfaces. It has been shown Ogura  that In exchanges at Al acid sites of the zeolite and, due to its specific surface, 90 % of these sites are inside the channels. In/ZSM5 is known as a promising catalyst in one of the major topics in environmental catalysis, the selective catalytic reduction (SCR) of nitric oxides by hydrocarbons Ogura ; Miro . The nature of the active sites for SCR with methane (SCRM) in In/zeolite catalyst has been extensively studied but only little qualitative information about them has been obtainedArmor . Ogura et al. Ogura  showed that the active site in In/ZSM5 for the SCRM are intrapore In species, and suggested that they are (InO)+ oxoions coordinated at the Z- exchange sites of the zeolite. Very recently, Schmidt et al.Schmidt  reported an exhaustive spectroscopic study of the structure of In-species (not necessarily active) using several preparation methods to synthesize In/ZSM5. Since they did not succeed to obtain a sample with an isolated active species they have not attempted to extract a particular In-O distance among several oxygen neighbors to verify the existence of the In-oxo site. In a previous work we studied in detail the several In species that can arise using different preparation methods and activation treatments, leading to highly active catalysts Miro . By means of the PAC technique we characterized this system and quantified the amount of In active sites. However, little information about the structure of this species could be obtained because of the high distribution of its hyperfine frequency and the lack of accurate EFG calculations. In this work we have therefore selected the preparation route that yields to the largest amount of the active species and with the rest of In atoms forming bulk In2O3 at the external surface of the zeolite.

NH4-ZSM5 with a Si/Al ratio of 26.4 obtained as in Ref.Miro,  was used as starting material. Indium was incorporated to the zeolite by the conventional wet impregnation method, stirring an aqueous solution of InCl3 (added in an amount as to obtain a sample with 4 wt$`\%`$ of In) and NH4-ZSM5 at 80oC until all water was evaporated, followed by drying in a stove at 120oC for 12 h. After this, the solid was pretreated in a dried oxygen atmosphere heating up to 500oC at 5oC/min, and holding the final temperature for 12 h. Afterwards, the sample was calcined for 2 h in oxygen at 750oC. In order to perform the PAC experiments the probe 111In was introduced by adding traces of 111InCl3 to the non-radioactive InCl3 solution. 111In decays by electron-capture to excited states of 111Cd, which reaches the ground state mainly through a $`\gamma `$-$`\gamma `$ cascade with an intermediate nuclear state of quadrupole moment Q that interacts with the EFG existing at the probe site. This EFG, which perturbs the $`\gamma `$-$`\gamma `$ angular correlation, can be determined by measuring the time modulation of the $`\gamma `$-$`\gamma `$ coincidence rate. A four BaF2-detector fast-fast coincidence system in a coplanar 90o arrangement was used to measure the $`\gamma `$-$`\gamma `$ coincidence spectra, which are combined to form the anisotropy ratio R(t). This R(t) can be expressed as the sum of cosine functions of three interaction frequencies from which the strength, VZZ, and symmetry, $`\eta =(V_{XX}V_{YY})/V_{ZZ}`$ , of the diagonalized and traceless EFG tensor can be determined. If the probes are located at different sites, their relative concentration fi can also be determined model independently. A detailed description of the technique can be found elsewhere Frauenfelder .

Figure 1 shows the PAC spectrum for the In(111In)/ZSM5 sample taken at 500oC in air, and its corresponding Fourier transform. It should be noted that no after-effects due to the electron-capture decay of 111In are present, as expected at this temperature Habenicht . In Table I, the fitted hyperfine parameters f, $`\omega _Q=eQV_{ZZ}/40\mathrm{}`$, $`\eta `$ and $`\delta `$ (that accounts for a Lorentzian EFG distribution around a mean $`\omega _Q`$ value originated from a distribution of very similar neighborhoods of the probe for a certain site) are shown. The three hyperfine interactions indicate that In occupies three different sites. The first two (I1 and I2) correspond to 111In in the two inequivalent sites of In2O3 present in the expected 3:1 population ratio Habenicht , and the third one (I3) with a concentration of 46% corresponds to the catalytically active species, since it has been shown that In2O3 is not active for the reaction of interest Miro . The low $`\delta `$ values of I1 and I2 agree well with the fact that crystallites of very small dimension are really a minority in the sample. The slight deviation from axial symmetry of site D has been already reported for crystalline In2O3Bibiloni .

The In K-edge (27940 eV) X-ray absorption coefficients were measured at the D42 beamline (XAS13 station) of the DCI synchrotron at LURE. The spectra were recorded with 4 eV steps at room temperature in air using a two-crystals Ge400 monochromator and Kr gas in the ionization chambers. The EXAFS data were extracted from the measured absorption spectra by standard methods Cook . The spectra were Fourier-transformed in order to obtain the radial distribution functions around In atoms in the samples. This function gives a view of the atomic distribution around the absorber, with peaks at distances where neighboring shells are located. The height of the peaks is related to the type and coordination number (CN) of neighbors. The study of small particles by EXAFS could lead to average coordination numbers (ACN) smaller than those expected for the bulk because of the contribution of atoms in the surface. The structural and thermal disorder, and a distribution of slightly different surroundings for a certain site (accounted for by the Debye-Waller factor, $`\sigma ^2`$) could also reduce the height of the peaks.

Figure 2 shows the K-edge EXAFS spectra of In2O3 and In/ZSM5 as well as their Fourier transforms. At first sight, both spectra differ only in the smaller amplitude of the second one. As there is a mixture of indium species in the catalyst - In2O3 and In-oxo active species - and both of them have In-O bonds with overlapping bond-length, it is not possible to fit directly the EXAFS spectra to obtain a deconvoluted information of one of them. In contrast, if doing so we would obtain an average of the In CN of both species weighted by the relative abundance of each species in the catalyst. Even though it should be possible to propose a two-shells fitting of the spectrum with modified CN’s, this procedure will not lead to a unique solution. At this point, the relative fractions of each species obtained from 111In-PAC measurements in a model-independent way could be used to isolate the active species by subtracting the spectrum of In2O3 weighted with its relative fraction, a method similar to the difference technique used to isolate a minor component from the major components in an EXAFS spectrum Teo . Nevertheless, there is one thing that must be considered, i.e., the size of the In2O3 crystallites. The $`\eta `$ parameter of the symmetric site D observed in the PAC experiments suggested that small crystallites of In2O3 could be present.

Figure 3 shows the ACN of In at the symmetric site of neutral spherical particles of In2O3 for different particle sizes, obtained averaging the CN of each In atom in the particle. As can be seen, the ACN grows rapidly for small particles and approaches the CN of the bulk (CN=6) for radii greater than ca. 20 Å. This implies that, if the In2O3 particles present in our sample are very small (with radii of the order of 5 to 10 Å), the subtraction of the In2O3 bulk spectrum to the In/ZSM5 one would lead to an overestimation of the In CN and hence to an underestimation of the CN of the active species. On the contrary, if the sesquioxide particles are bigger than 20 Å in radius, their ACN would be larger than 5.5 and it could be correct, in a first approximation, to subtract the In2O3 bulk spectrum to obtain that of the other species present in the catalyst.

Hence, the SAXS technique was used to determine the volume distribution of indium oxide particles in the In/ZSM5 catalysts. SAXS experiments were performed at the SAS beamline Kellerman  at the National Synchrotron Light Laboratory (LNLS-Campinas, Brazil). The cluster volume distribution function DV(R) Glatter  of the In-Ox particles is shown in Fig.4.

This function was determined from the SAXS intensity profile (Fig. 4, inset) using the GNOM package Semenyuk . This profile was obtained as usual Benedetti  by subtracting the ZSM5 normalized SAXS spectrum to that of the In/ZSM5 catalysts to get rid of the scattering contribution coming from the porous ZSM5 support. Figure 4 shows that there is a broad distribution of In2O3 particle sizes in the catalysts and there is a maximum of DV(R) at R=50 Å with a long tail to larger radii. Although there are small particles with R $`<`$ 20 Å, its relative fraction is very small. Partial integration of the curve shows that 93% of the volume corresponds to particles with radius bigger than 20 Å.

For particles bigger than 20 Å in radius, the average coordination number is close to that of the bulk and the difference can be included within the experimental error for the CN (ca. 10%). Under these circumstances, it is valid to subtract the EXAFS spectrum of bulk In2O3 to the total signal obtained for the catalyst to obtain the EXAFS spectrum of the I3 site.

Figure 5 shows the Fourier transform of the difference spectrum obtained subtracting the In2O3 spectrum weighted with the percentage (54%) of the oxide in the catalysts - determined by PAC - to the In/ZSM5 EXAFS spectrum. The inset shows the backtransform of the peak in the range 0.9-2.6 Å and the fitted EXAFS functions obtained using phase and amplitudes generated by the FEFF code Zabinsky . The parameters obtained are shown in Table II. Two oxygen shells and one aluminium shell were found in the fitting procedure. One oxygen atom is located at 2.10 Å while 2 oxygen atoms are at 2.20 Å. These last two oxygen atoms would belong to the zeolite structure and would be bound to an aluminum atom forming an acid site. The distance of the aluminum atom to the In atom is 2.59 Å. The closer oxygen atom is forming the (InO)+ species proposed as the active site and the shorter distance is consistent with a double bond. It is seen that the Debye-Waller factor of the second and third shells are large. This behavior can be understood because of the different configurations in which the (InO)+ species can be bound to different acid sites of the zeolite, leading to slightly different distances In-O (2nd shell) and In-Al. Similar $`\sigma `$ values have been found for small Pt clusters inside zeolite channels Jentys . Using theoretical calculations, Jentys et al. found that Pt particles inside zeolite channels have Debye-Waller factors two orders of magnitude larger than an isolated Pt particle. They assigned this effect to an exponential damping of the EXAFS signal resulting from contributions of the zeolite lattice atoms. The several anchoring configurations of the (InO)+ species would also explain the high EFG distribution of I3 found by PAC (see Table I). Because of the size of the monometallic species, the active site can only be found in the 10-atom ring (8 Å diam.) of the ZSM5 zeolite as it would not fit in the 6-atom ring (5 Å diam.).

Consequently, by means of XAFS, TDPAC and SAXS experiments we determined the local structure of extremely nanozised In-O species (with the geometry indicated in Table II), responsible for the SCR of NOx, located inside the larger ZSM5 zeolite channel.

We wish to dedicate this article to our colleague and friend Patricia Massolo in the 10th anniversary of her early death. This research was partially supported by ANPCyT, CONICET, and Fundación Antorchas (Argentina), by LURE (CNRS, France), and LNLS (Brasil) under project SAXS 314/97. Authors express gratitude to I. Torriani for her assistance at the SAXS beamline.