text:
# The nuclear quadrupole moment of 57Fe from microscopic nuclear and atomic calculations

\[

## Abstract

The nuclear quadrupole moment of the $`I^\pi `$=3/2- excited nuclear state of 57Fe at 14.41 keV, important in Mössbauer spectroscopy, is determined from the large-scale nuclear shell-model calculations for 57Fe and also from the electronic ab initio and density functional theory calculations including solid state and electron correlation effects for the molecules Fe(CO)5 and Fe(C5H5)2. Both independent methods yield very similar results. The recommended value is 0.16(1) eb. The NQM of the isomeric 10+ in 54Fe has also been calculated. The new value (0.5 eb), consistent with the perturbed angular distribution data, is by a factor of two larger than the currently recommended value.

\]

Mössbauer spectroscopy of 57Fe plays an important role in the structural determination of iron containing solid state compounds. In principle, the nuclear quadrupole moment (NQM) of the isomeric $`I`$=3/2 state in 57Fe can be determined from Mössbauer data; however, the analysis requires the calculation of the electric field gradient (EFG). As these atomic calculations are quite involved, studies employing different methods arrived at quite distinct results; the values of NQM in the range from –0.19 to +0.44 eb have been reported . This quite unsatisfying situation could also not have been settled by nuclear structure calculations of the NQM as calculations within the nuclear shell model, the most reliable tool for such studies, had to be performed in strongly truncated model spaces and with rather untested effective interactions. In recent years, decisive progress has been achieved in both the atomic calculations of the EFG and in nuclear shell model studies.

In 1995, Dufek and co-workers applied the density functional theory for a series of iron-containing solid state compounds. For 57Fe they obtained a NQM of 0.16 eb , in contradiction with the previously accepted value of 0.082 eb obtained from Hartree-Fock (HF) EFG calculations  and from truncated nuclear shell-model calculations combined with the perturbed angular distribution data . In subsequent work, Su and Coppens obtained a NQM value of 0.12(3) eb using Sternheimer-corrected EFGs . In this Letter, we shall demonstrate that state-of-the-art nuclear and atomic physics calculations lead to the same NQM for the 57Fe isomeric state, settling a long-standing controversial issue.

We shall begin with the nuclear physics discussion. The nucleus 57Fe has two low-lying $`3/2^{}`$ states which are experimentally split by only 353 keV. To describe the structure of such nearly degenerate states is quite demanding. Clearly, the large-scale shell model is the method of choice. Due to recent progress in programming and hardware development, modern shell-model calculations based on microscopic effective interactions can handle configuration spaces that were prohibitively large only several years ago . More specifically, modern diagonalization shell-model codes can now handle medium-mass nuclei ($`A`$=50-60) in the middle of the $`pf`$ shell in full $`0\mathrm{}\omega `$ space. To put things in perspective, the shell-model calculations of the Utrecht Group , which were used in Ref.  to extract the NQM of 57Fe, restricted the number of holes in the $`f_{7/2}`$ orbit to three.

Shell-model calculations depend crucially on two factors: the model space and the effective interaction. Our calculations have been performed using the code NATHAN . NATHAN has been developed in the $`jj`$-coupling scheme using quasispin formalism. We adopt a version of the code adapted to shared-memory parallel machines. For the two lowest $`3/2^{}`$ levels, we assumed the truncated space (containing 8,120,105 $`I^\pi `$=3/2- states) in which maximally 6 nucleons were allowed to be excited from the $`f_{7/2}`$ orbital to the rest of the $`pf`$ shell. To test the convergence of our results, we have performed a full pf shell calculation for the lowest $`3/2^{}`$ state. This calculation includes 25,743,302 states ($``$2$``$1013 non-zero matrix elements) and is one of the largest shell-model diagonalization performed to date. For all adopted interactions, the results of the truncated and complete calculations are identical.

During the last few years, a considerable effort went into the development of effective interactions in the $`pf`$ shell. In this work three different effective interactions have been employed: KB3F, KB3G, and FDP6. The interactions KB3F  and KB3G  are both reasonable attempts to correct the defects of the well-known KB3 interaction  in the upper part of the $`pf`$-shell. They differ slightly in their collectivity of states around 58Ni. As far as 57Fe is concerned, KB3G is clearly the interaction of choice; it has been shown to be very successful in describing experimental data (including energy levels and electromagnetic properties) in the mass region $`A`$=50-52 and around $`A`$=56 . The FDP6 interaction  was originally fitted to the spectroscopic properties of $`f_{7/2}`$ nuclei and has since been extended to nuclei around the $`N`$=$`Z`$=28 shell closure . For the effective charges, we took the quadrupole charges 1.5$`e`$ for protons and 0.5$`e`$ for neutrons when calculating E2 transitions and moments, and the spin and orbital gyromagnetic factors $`g^s`$=0.75$`g_{\mathrm{bare}}^s`$, $`g_\pi ^l`$=1.1 $`\mu _N`$, and $`g_\nu ^l`$=$``$0.1 $`\mu _N`$ for M1 transitions and moments, e.g. . We adopt the oscillator parameter $`b=1.01A^{1/6}`$ fm.

Our shell-model calculations reproduce well the experimental level scheme of 57Fe. In particular, the low-lying $`I^\pi `$=3/2- doublet is predicted by theory. Experimentally, the energy difference between these states is 353 keV, while theoretically it is 526, 209, and 465 keV using the KB3F, KB3G, and FPD6 interactions, respectively. The third 3/2- state appears much higher in energy, at about 1.6 MeV. Table I summarizes our shell-model results for the quadrupole and magnetic moments of the two lowest $`3/2^{}`$ states in 57Fe. It is seen that these close-lying states in the doublet are predicted to have very different shapes and magnetic moments, i.e., their E2 and M1 moments have different signs. Consequently, even small changes in the shell-model interaction, hence in the coupling between these states, significantly impact theoretical predictions. We note that a convergence of the shell-model results is not reached if only three $`f_{7/2}`$ holes are allowed, as was done in the previous truncated calculations . For example, the value of the NQM of the 3/2$`{}_{1}{}^{}{}_{}{}^{}`$ state, $`Q(3/2_1^{})`$, obtained in such truncated calculations with the KB3G interaction, is 0.097 eb, as compared to the exact value of 0.064 eb.

In order to assess the sensitivity of the results on the adopted effective force, we computed overlaps between wave functions of the $`3/2^{}`$ states obtained with the different interactions. It turned out that the spaces of the two lowest $`3/2^{}`$ states in KB3F and KB3G are practically (up to 95%) the same. That is, both $`3/2^{}`$ wave functions obtained with KB3F can be approximately derived from the $`3/2^{}`$ doublet calculated with KB3G by simple rotation. This does not hold for KB3G and FDP6. Since FPD6 tends to excite particles to the $`f_{5/2}`$ shell rather than the $`p_{3/2}`$ shell (FPD6 puts the $`f_{5/2}`$ orbit too low at $`N`$=28), the contributions from higher-lying $`3/2^{}`$ states amount to about 37%.

Due to the great sensitivity of the shell-model predictions caused by the not-very-well-controlled off-diagonal coupling between the $`3/2^{}`$ doublet, the calculations need to be constrained by the available experimental data. To this end, we choose the magnetic moment of the $`3/2_1^{}`$ state, which has been precisely determined experimentally, $`\mu (3/2_1^{})`$=$``$0.1549(2) $`\mu _N`$ Assuming that the “true” wave functions of the $`3/2^{}`$ doublet are given by a simple rotation of the shell-model states,

$$|3/2^{}=\alpha _1|3/2_1^{};\mathrm{SM}+\alpha _2|3/2_2^{};\mathrm{SM},$$
(1)
one can determine the mixing amplitudes $`\alpha _1`$ and $`\alpha _2`$ ($`\alpha _1^2+\alpha _2^2=1`$) by requesting that the measured value of $`\mu (3/2_1^{})`$ be reproduced. The results of such two-level mixing calculations are displayed in Table II. Contrary to the pure shell-model results, one obtains a satisfying agreement between KB3F and KB3G, and the data. In particular, the quadrupole moments predicted by these interactions are very similar, $`Q(3/2_1^{})`$$``$0.17 eb. The only serious difference is the sign of the magnetic moment of the $`3/2_2^{}`$ state. Unfortunately, this quantity is poorly determined experimentally: $`\mu (3/2_2^{})`$$`<`$0.6 $`\mu _N`$ The predictions of FPD6 for the E2 moments are rather far from the data: both $`3/2^{}`$ states are calculated to be practically spherical, and the $`B(E2)`$ transition connecting these states is enhanced by a factor of $``$40. Part of this failure comes from too-low a placement of the $`f_{5/2}`$ level in FPD6, as discussed in detail in Ref. , where more examples of KB3G and FPD6 calculations can be found. Therefore, we conclude that the best shell-model estimate of the quadrupole moment of the $`3/2_1^{}`$ state in 57Fe is 0.16 eb.

In Ref.  it has been pointed out that the ratio of quadrupole moments of the isomeric $`10^+`$ state in 54Fe to the $`3/2^{}`$ state in 57Fe, $`Q[^{54}`$Fe$`(10^+)]`$/$`Q[^{57}`$Fe$`(3/2_1^{})]`$, represents a stringent constraint. Possible uncertainties in the atomic EFG calculations cancel largely in the experimental determination of the ratio from Mössbauer data. The most recent value is $`Q[^{54}`$Fe$`(10^+)]`$/$`Q[^{57}`$Fe$`(3/2_1^{})]`$ =3.62$`\pm `$ 0.22 . Unfortunately, the NQM of the isomeric state in 54Fe is not known experimentally. We have thus calculated the NQM of the isomeric $`10^+`$ state in 54Fe in a large-scale shell-model study in which we allowed 8 particles to be promoted from the $`f_{7/2}`$ orbital to the rest of the shell. We have checked the convergence of our results by performing calculations in which only 6 particles could be promoted; the NQMs were the same at both levels of truncation. For all the interactions used, we find rather similar NQMs: 0.51 eb (KB3F), 0.50 eb (KB3G), and 0.56 eb (FDP6). Using the NQMs of Table II for the isomeric state in 57Fe, these values correspond to $`Q[^{54}`$Fe$`(10^+)]`$/$`Q[^{57}`$Fe$`(3/2_1^{})]`$ ratios which agree rather well with the measurement for the KB3F and KB3G interactions (2.9 and 3.1, respectively). For the FDP6 interaction, the ratio is significantly too large, as the NQM of the isomer in 57Fe is severely underestimated. We note that our predicted values for the magnetic moment of the $`10^+`$ state in 54Fe (6.4 $`\mu _N`$ in KB3F, 6.5 $`\mu _N`$ in KB3G, and 7.1 $`\mu _N`$ in FDP6) compare well with the experimental value of 7.28$`\pm 0.01`$ n.m .

We now turn to our atomic physics studies. The use of the density functional theory for the calculation of EFGs in transition metals is questionable. Recent calculations on CuF showed a variation of the Cu EFG ranging from +0.50 au for the local density approximation (LDA), and +0.44 au at the generalized gradient approximation (GGA) level (BPW91) to +0.07 au at the hybrid level (B3PW91) , as compared to the experimental value of $``$0.31(2) au . In contrast, relativistic ab initio coupled-cluster calculations give $``$0.34 au , in perfect agreement with the experimental result.

In this work, we carried out density functional as well as ab initio calculations for the molecules Fe(CO)5 and Fe(C5H$`{}_{5}{}^{})_2`$. We have adopted a wide range of exchange and correlation functionals for the electronic structure calculations of the free molecules Fe(CO)5 and Fe(C5H$`{}_{5}{}^{})_2`$ (for the terminology see Refs. ): Hartree-Fock-Slater (HFS), X$`\alpha `$, Local Density Approximation (LDA), the GGA functionals B-HFS (Slater exchange plus Becke nonlocal exchange), B-LYP (B-HFS plus the correlation functional of Lee, Yang and Parr), B-PW91 (B-HFS plus the correlation functional of Perdew and Wang), the hybrid GGA functionals B3-LYP (Becke three-parameter functional), B3-PW91 (same as B3-LYP except with the nonlocal correlation term of Perdew and Wang), BHH (Becke half-and-half together with the LYP correlation functional), and BHH-LYP (same as BHH but with 0.5 of the Becke nonlocal exchange term added to the energy). For comparison with DFT, we carried out the HF many-body perturbation theory (MBPTn) up to order n=4, as well as coupled-cluster singles doubles including non-iterative triple excitations \[CCSD(T)\] . The electronic coupled cluster calculations required three months of CPU time and 20 Gbytes of disk storage on a Origin 2000 SGI. For Fe(CO)5 we investigated solid state effects to the iron electric field gradient by performing HF, B3LYP, and LDA calculations using the program CRYSTAL98  and the solid state structure given by Böse and Bläser . Detailed structural data and basis sets used will be published elsewhere .

The results of our electronic structure calculations for the iron EFGs of Fe(CO)5 and Fe(C5H$`{}_{5}{}^{})_2`$ are shown in Table III. We first note that the single-reference many-body perturbation theory shows extreme oscillatory behavior and is practically useless for the determination of transition element EFGs. DFT results range from 1.09 au to 1.57 au for Fe(CO)5, and from 1.36 au to 2.49 au for Fe(C5H$`{}_{5}{}^{})_2`$, depending on the functional applied. If we accept the coupled-cluster result as the most accurate value, the hybrid GGA functionals perform well for Fe(CO)5, while the non-hybrid GGA functionals are preferred for Fe(C5H$`{}_{5}{}^{})_2`$. We note that for CuCl the BHH functionals gave the best description . This is clearly not a satisfying situation. On the other hand, an encouraging result is that solid state effects from nearest-neighbor interactions in Fe(CO)5 are small and can basically be neglected. Based on the CCSD(T) EFGs, we obtain from the Mössbauer data of Fe(CO)5) (2.57 mm/sec ) and Fe(C5H$`{}_{5}{}^{})_2`$ (2.4 mm/sec ) a NQM of 0.177 eb and 0.159 eb, respectively. However, basis set incompleteness and relativistic effects may increase the iron EFGs and, therefore, further decrease the NQM . Consequently, our best estimate using EFGs together with Mössbauer data for these molecules is 0.15(2) eb. In order to improve further on these results, large-scale relativistic coupled-cluster calculations are necessary, which are currently not feasible for such big molecules. We emphasize, however, that EFGs obtained from current DFT for transition element compounds should be taken with some care as the results in Table III show.

In summary, the quadrupole moment of the first $`3/2^{}`$ state in 57Fe at 14.41 keV, important in Mössbauer spectroscopy, has been calculated using state-of-the-art nuclear and atomic models. Both calculations yield results which are consistent with $`Q(3/2_1^{})`$$``$0.16 eb, in nice agreement with the recently reported value  and the results of Ref. . As a by-product of this work, the NQM of the isomeric 10+ in 54Fe has also been calculated. The new value, consistent with the perturbed angular distribution data, is 0.5 eb, i.e., a factor of two greater than the currently recommended value .

This research was supported by the U.S. Department of Energy under Contract Nos. DE-FG02-96ER40963 (University of Tennessee), DE-FG02-97ER41019 (University of North Carolina), DE-AC05-00OR22725 with UT-Battelle, LLC (Oak Ridge National Laboratory), DE-FG02-92ER40694 (Tennessee Technological University), the Danish Research Foundation, the Schweizer Nationalfond (grant 20-61822.00), the Marsden Fund Wellington, the Royal Society of New Zealand, and the Auckland University Research Committee. The Auckland group is indebted to the Alexander von Humboldt Foundation for financial support. GMP wants to thank the Carlsberg Foundation for a fellowship.


format: text

text:
# Bessel Integrals and Fundamental Solutions for a Generalized Tricomi Operator

## 1 Introduction

Consider the generalized Tricomi operator

(1.1)
$$P=y\mathrm{\Delta }+\frac{^2}{y^2},$$

where $`\mathrm{\Delta }=_{j=1}^n{\displaystyle \frac{^2}{x_j^2}}`$ is the Laplace operator. Our aim is to find, via partial Fourier transform with respect to $`x=(x_1,\mathrm{},x_n)`$, fundamental solutions relative to an arbitrary point $`(a,0)`$ on the hyperplane $`y=0`$ in $`^{n+1},`$ that is, distributions that are solutions to

$$Pu=\delta (xa,y),$$

where $`\delta (xa,y)`$ is the Dirac measure concentrated at $`(a,0).`$ Since $`P`$ is invariant under translations parallel to that hyperplane, it suffices to consider the case when the Dirac measure is concentrated at the origin.

If $`n=1,`$ then (1.1) is the classical Tricomi operator,

(1.2)
$$𝒯=y\frac{^2}{x^2}+\frac{^2}{y^2}.$$

It is known that, for this operator, the equation $`9x^2+4y^3=0`$ defines the two characteristic curves that originate at the origin. They divide the plane in two disjoint regions, namely,

$$D_+=\{(x,y)^2:9x^2+4y^3>0\},$$

the region “outside” the characteristics and

$$D_{}=\{(x,y)^2:9x^2+4y^3<0\},$$

the region “inside” the characteristics. Note that $`D_{}`$ is entirely contained in the hyperbolic region $`y<0.`$

In the paper  it was shown the existence of the following two fundamental solutions relative to the origin:

(1.3)
$$F_+(x,y)=\{\begin{array}{cc}\frac{\mathrm{\Gamma }(1/6)}{32^{2/3}\pi ^{1/2}\mathrm{\Gamma }(2/3)}(9x^2+4y^3)^{1/6}\hfill & \text{in }D_+\hfill \\ & \\ 0\hfill & \text{elsewhere}\hfill \end{array}$$

and

(1.4)
$$F_{}(x,y)=\{\begin{array}{cc}\frac{3\mathrm{\Gamma }(4/3)}{2^{2/3}\pi ^{1/2}\mathrm{\Gamma }(5/6)}|9x^2+4y^3|^{1/6}\hfill & \text{in }D_{}\hfill \\ & \\ 0\hfill & \text{elsewhere}.\hfill \end{array}$$

The support of $`F_{},`$ the closure of $`D_{},`$ is, except for the origin, entirely contained in the hyperbolic region ($`y<0`$), while the support of $`F_+,`$ the closure of $`D_+,`$ consists of the whole elliptic region ($`y>0`$), the parabolic region ($`x`$-axis), and extends to the hyperbolic region up to and included the characteristic curves. The method used in  to prove these results was based upon the property that solutions to the equation $`𝒯u=0`$ are invariant under the dilation $`d_t(x,y)=(t^3x,t^2y)`$ in $`^2.`$

In  we proved formulas (1.3) and (1.4) by using partial Fourier transform with respect to the $`x`$ variable. In this paper we extend the results of  to the generalized Tricomi operator (1.1). Since the dimension of the space variable $`x`$ is now $`n>1,`$ technical difficulties in evaluating Fourier transforms involving Bessel functions do occur. We circumvent them by calculating integrals of the type

$$I_ϵ(a,b)=_0^{\mathrm{}}e^{ϵt}t^\lambda J_\mu (at)J_\nu (bt)𝑑t$$

(Section 3). As a result, we obtain for the operator (1.1) the following fundamental solutions:

$$F_+(x,y)=\{\begin{array}{cc}\frac{3^{n2}}{2^{2/3}\pi ^{n/2}}\frac{\mathrm{\Gamma }(\frac{n}{2}\frac{1}{3})}{\mathrm{\Gamma }(2/3)}(9|x|^2+4y^3)^{\frac{1}{3}\frac{n}{2}}\hfill & \text{in }D_+^n\hfill \\ & \\ 0\hfill & \text{elsewhere},\hfill \end{array}$$

whose support is the closure of the region $`D_+^n=\{(x,y)^{n+1}:9|x|^2+4y^3>0\}`$ and

$$F_{}(x,y)=\{\begin{array}{cc}\frac{3^n\mathrm{\Gamma }(\frac{4}{3})}{2^{2/3}\pi ^{n/2}\mathrm{\Gamma }(\frac{4}{3}\frac{n}{2})}|9|x|^2+4y^3|^{\frac{1}{3}\frac{n}{2}}\hfill & \text{in }D_{}^n\hfill \\ & \\ 0\hfill & \text{elsewhere}\hfill \end{array}$$

supported by the closure of the region $`D_{}^n=\{(x,y)^{n+1}:9|x|^2+4y^3<0\}.`$

These fundamental solutions clearly generalize formulas (1.3) and (1.4). In addition, we also show that no matter the parity of the dimension $`n,`$ $`F_{}(x,y)`$ does not have support on the boundary of the region $`D_{}`$ in the hyperbolic region and hence, in contrast with what happens for strictly hyperbolic operators – like the wave operator – the Huyghens principle does not hold for the generalized Tricomi operator.

It is well known  that, for the wave operator, the regularity of the fundamental solutions degenerates as the dimension $`n`$ increases: a locally constant function, when $`n=1;`$ an absolutely continuous measure relative to the Lebesgue measure, when $`n=2;`$ a measure carried by the surface of the forward light-cone, when $`n=3;`$ and so on. However, for the generalized Tricomi operator they always remain locally integrable.

The plan of this paper is as follows. In Section 2 we reduce, via partial Fourier transform, the original problem to an equivalent one of finding fundamental solutions for a second order ordinary differential equation and show how these can be represented in terms of Airy functions or Bessel functions. In Section 3 we obtain explicit formulas for Fourier transforms of the functions $`|\xi |^{\pm \nu }J_\nu (|\xi |),`$ $`|\xi |^\nu K_\nu (|\xi |),`$ and $`|\xi |^\nu N_\nu (|\xi |),`$ $`\xi =(\xi _1,\mathrm{},\xi _n),`$ $`n1.`$ These formulas are used to construct the fundamental solution supported by $`\overline{D_{}^n}`$ (Section 4) and the fundamental solution supported by $`\overline{D_+^n}`$ (Section 5). In the Appendix the reader will find the definitions of the Bessel functions $`J_\nu (z),`$ $`I_\nu (z),`$ $`K_\nu (z),`$ and $`N_\nu (z),`$ and a list of properties that are needed throughout this work. We also recall the definition of hypergeometric functions and some of their main properties.

The method of finding solutions, via partial Fourier transform, to the equation

(1.5)
$$Pu=y\mathrm{\Delta }u+\frac{^2u}{y^2}=f$$

with, say, $`f𝒞_c^{\mathrm{}}(^{n+1})`$ is a natural one although not particularly new. In the monograph , R. J. P. Groothuizen exhibits a solution $`u`$ in terms of a Fourier integral operator whose symbol is obtained by the partial Fourier transform analysis employed in this paper. He also obtains fundamental solutions given by Fourier integral operators. Our results are more precise since we give explicit formulas for the fundamental solutions and these are tempered distributions. Consequently, we can obtain more transparent representations for the solution $`u`$ to the equation (1.5) as a convolution of $`f`$ with any of the fundamental solutions here described.

## 2 Preliminaries

Consider the more general problem of finding fundamental solutions for the operator (1.1) relative to an arbitrary point $`(0,b)`$ on the $`y`$-axis. That is, one wishes to find distributional solutions to the equation

(2.1)
$$y\mathrm{\Delta }F+\frac{^2}{y^2}F=\delta (x)\delta (yb).$$

Partial Fourier transform with respect to $`x`$ reduces this problem into finding fundamental solutions to the ordinary differential equation

(2.2)
$$\stackrel{~}{F}_{yy}y|\xi |^2\stackrel{~}{F}=\delta (yb),$$

where $`\stackrel{~}{F}`$ denotes the partial Fourier transform of $`F,`$ $`|\xi |^2=_{j=1}^n\xi _j^2,`$ and $`\delta (yb)`$ is the Dirac measure concentrated at $`b.`$

One way of solving (2.2) (see ) consists in selecting two linearly independent solutions $`U_1(\xi ,y)`$ and $`U_2(\xi ,y)`$ to the homogeneous equation

(2.3)
$$\stackrel{~}{u}_{yy}y|\xi |^2\stackrel{~}{u}=0$$

so that their Wronskian at $`y=b`$ (which in the case under consideration is the same as the Wronskian at $`y=0`$) is equal to $`1`$ and in defining

(2.4)
$$\stackrel{~}{F}(\xi ,y;b)=\{\begin{array}{cc}U_2(\xi ,b)U_1(\xi ,y)\hfill & \text{if }yb\hfill \\ & \\ U_1(\xi ,b)U_2(\xi ,y)\hfill & \text{if }yb.\hfill \end{array}$$

It is a matter of verification that $`\stackrel{~}{F}(\xi ,y;b)`$ is a solution to (2.2).

Equivalently , one can obtain a fundamental solution to the equation (2.2) by finding two linearly independent solutions to the homogeneous equation (2.3): $`u_1(\xi ,y;b),`$ defined for $`y>b,`$ and $`u_2(\xi ,y;b),`$ defined for $`y<b,`$ so that

(i) the limit of $`u_1`$ as $`yb+`$ equals the limit of $`u_2`$ as $`yb`$

and

(ii) the limit of $`u_{1,y}`$ as $`yb+`$ minus the limit of $`u_{2,y}`$ as $`yb`$ is equal to 1. The function $`\stackrel{~}{F}(\xi ,y;b)`$ is now defined by

(2.5)
$$\stackrel{~}{F}(\xi ,y;b)=\{\begin{array}{cc}u_1(\xi ,y;b)\hfill & \text{if }yb\hfill \\ & \\ u_2(\xi ,y;b)\hfill & \text{if }yb.\hfill \end{array}$$

In what follows we will use at our convenience either one of these two expresions for $`\stackrel{~}{F}(\xi ,y;b).`$

Returning to the original problem (2.1), if we can choose $`U_1(\xi ,y)`$ and $`U_2(\xi ,y),`$ in formula (2.4), or $`u_1(\xi ,y;b)`$ and $`u_2(\xi ,y;b),`$ in formula (2.5), so that $`\stackrel{~}{F}(\xi ,y;b)`$ is a tempered distribution with respect to $`\xi =(\xi _1,\mathrm{},\xi _n),`$ then its inverse Fourier transform of $`\stackrel{~}{F}(\xi ,y;b)`$ defines a fundamental solution to the operator (1.1).

We now proceed to find linearly independent solutions to the homogeneous equation (2.3). The change of variables $`z=|\xi |^{2/3}y`$ transforms (2.3) into the classical Airy’s equation $`u^{\prime \prime }zu=0.`$ Two linearly independent solutions to that equation are $`Ai(z)`$ and $`Bi(z)`$ respectively called the Airy functions of the first and second kinds. It is known  that these two functions can be represented in terms of Bessel functions of order $`\pm 1/3`$ as follows. If $`|\mathrm{arg}z|<(2\pi /3),`$ then

(2.6)
$$\{\begin{array}{cc}Ai(z)=\frac{z^{1/2}}{3}[I_{1/3}(\frac{2}{3}z^{3/2})I_{1/3}(\frac{2}{3}z^{3/2})]=\frac{1}{\pi }(\frac{z}{3})^{1/2}K_{1/3}(\frac{2}{3}z^{3/2}),\hfill & \\ & \\ Bi(z)=(\frac{z}{3})^{1/2}[I_{1/3}(\frac{2}{3}z^{3/2})+I_{1/3}(\frac{2}{3}z^{3/2})].\hfill & \end{array}$$

If $`|\mathrm{arg}z|<(2\pi /3),`$ then

(2.7)
$$\{\begin{array}{cc}Ai(z)=\frac{z^{1/2}}{3}[J_{1/3}(\frac{2}{3}z^{3/2})+J_{1/3}(\frac{2}{3}z^{3/2})],\hfill & \\ & \\ Bi(z)=(\frac{z}{3})^{1/2}[J_{1/3}(\frac{2}{3}z^{3/2})J_{1/3}(\frac{2}{3}z^{3/2})].\hfill & \end{array}$$

In the Appendix I, where a brief review of Bessel functions is presented, we show that the following relations hold:

(2.8)
$$Ai(0)=\frac{3^{2/3}}{\mathrm{\Gamma }(2/3)},Ai^{}(0)=\frac{3^{4/3}}{\mathrm{\Gamma }(4/3)}$$

and

(2.9)
$$Bi(0)=\frac{3^{1/6}}{\mathrm{\Gamma }(2/3)},Bi^{}(0)=\frac{3^{5/6}}{\mathrm{\Gamma }(4/3)}.$$

As a consequence, the Wronskian of $`Ai(z)`$ and $`Bi(z)`$ evaluated at $`z=0`$ is

(2.10)
$$W(Ai(z),Bi(z))_{|_{z=0}}=1/\pi .$$

Indeed, we have

$$W(Ai(z),Bi(z))_{|_{z=0}}=\left|\begin{array}{cc}\frac{3^{2/3}}{\mathrm{\Gamma }(2/3)}& \frac{3^{1/6}}{\mathrm{\Gamma }(2/3)}\\ & \\ \frac{3^{4/3}}{\mathrm{\Gamma }(4/3)}& \frac{3^{5/6}}{\mathrm{\Gamma }(4/3)}\end{array}\right|=\frac{23^{3/2}}{\mathrm{\Gamma }(2/3)\mathrm{\Gamma }(4/3)}=\frac{1}{\pi },$$

because

(2.11)
$$\mathrm{\Gamma }(2/3)\mathrm{\Gamma }(4/3)=\frac{2\pi }{3^{3/2}}.$$

We now choose the following two linearly independent solutions to the equation (2.3):

$$U_1(\xi ,y)=\sqrt{\pi }|\xi |^{1/3}Ai(|\xi |^{2/3}y)\text{and}U_2(\xi ,y)=\sqrt{\pi }|\xi |^{1/3}Bi(|\xi |^{2/3}y),$$

and note that, by virtue of (2.10), the Wronskian of $`U_1(\xi ,y)`$ and $`U_2(\xi ,y)`$ is equal to $`1.`$ Next, according to (2.4), the distribution

(2.12)
$$\stackrel{~}{F}(\xi ,y;b)=\{\begin{array}{cc}\pi |\xi |^{2/3}Bi(|\xi |^{2/3}b)Ai(|\xi |^{2/3}y)\hfill & \text{if }yb\hfill \\ & \\ \pi |\xi |^{2/3}Ai(|\xi |^{2/3}b)Bi(|\xi |^{2/3}y)\hfill & \text{if }yb\hfill \end{array}$$

is a fundamental solution to the ordinary differential equation (2.2). A fundamental solution $`F(x,y;b)`$ to the generalized Tricomi operator (1.1) and relative to the point $`(0,b)`$ would then be the inverse Fourier transform of $`\stackrel{~}{F}(\xi ,y;b)`$ whenever that Fourier transform exists.

We do not know how to obtain an explicit formula (or formulas) for the inverse Fourier transform of $`\stackrel{~}{F}(\xi ,y;b)`$ when $`b0,`$ a problem that merits to be investigated. We conjecture that when $`b<0,`$ that is, the point $`(0,b)`$ is in the hyperbolic region, there exists two fundamental solutions that converge, as $`b0,`$ to the fundamental solutions $`F_+(x,y)`$ and $`F_{}(x,y)`$ defined, respectively, by the formulas (5.6) and (4.2). The two fundamental solutions described in the monograph  do not seem to satisfy these requirements.

However, when $`b=0,`$ we will show in the following sections how to obtain from formula (2.12), as well as formulas similar to it, explicit expressions for fundamental solution to (1.1).

The case $`n=1`$. For sake of completeness and in order to motivate our work in the forthcoming sections we briefly present the results of the paper , for the Tricomi operator (1.2), that is the case when $`n=1.`$ If $`b=0`$ and we take into account the values of $`Ai(0)`$ and $`Bi(0)`$ as given by (2.8) and (2.9), then formula (2.12) simplifies as follows:

(2.13)
$$\stackrel{~}{F}(\xi ,y)=\{\begin{array}{cc}\frac{\pi |\xi |^{2/3}}{3^{1/6}\mathrm{\Gamma }(2/3)}Ai(|\xi |^{2/3}y)\hfill & \text{if }y0\hfill \\ & \\ \frac{\pi |\xi |^{2/3}}{3^{2/3}\mathrm{\Gamma }(2/3)}Bi(|\xi |^{2/3}y)\hfill & \text{if }y0,\hfill \end{array}$$

where, for simplicity, we wrote $`\stackrel{~}{F}(\xi ,y)`$ for $`\stackrel{~}{F}(\xi ,y;0).`$ Now the inverse Fourier transform of both expressions on the right-hand side of (2.13) can be explicitly calculated. To see this, first consider the representations of $`Ai(z)`$ and $`Bi(z)`$ in terms of Bessel functions as given by formulas (2.6) and (2.7):

$$Ai(|\xi |^{2/3}y)=\frac{1}{\pi }\left(\frac{|\xi |^{2/3}y}{3}\right)^{1/2}K_{1/3}(\frac{2}{3}|\xi |y^{3/2})$$

and

$$Bi(|\xi |^{2/3}y)=\left(\frac{|\xi |^{2/3}(y)}{3}\right)^{1/2}[J_{1/3}(\frac{2}{3}|\xi |(y)^{3/2})J_{1/3}(\frac{2}{3}|\xi |(y)^{3/2})].$$

Next by introducing the change of variables $`s=(2/3)y^{3/2},`$ whenever $`y0,`$ and $`t=(2/3)(y)^{3/2},`$ whenever $`y0,`$ rewrite formula (2.13) as

(2.14)
$$\stackrel{~}{F}(\xi ,y)=\{\begin{array}{cc}\alpha (\frac{s}{|\xi |})^{1/3}K_{1/3}(s|\xi |)\hfill & \text{if }y0\hfill \\ & \\ \beta (\frac{t}{|\xi |})^{1/3}[J_{1/3}(t|\xi |)J_{1/3}(t|\xi |)]\hfill & \text{if }y0\hfill \end{array}$$

where $`\alpha `$ and $`\beta `$ are constants given by

(2.15)
$$\alpha =\frac{1}{2^{1/3}3^{1/3}\mathrm{\Gamma }(2/3)}\text{and}\beta =\frac{\pi }{2^{1/3}3^{5/6}\mathrm{\Gamma }(2/3)}.$$

Theorem 3.1 of  proves that $`F(x,y),`$ the inverse Fourier transform of $`\stackrel{~}{F}(\xi ,y),`$ is then

$$F=\frac{3}{2}F_+\frac{1}{2}F_{},$$

a linear combination of the two fundamental solutions $`F_+`$ and $`F_{}`$ respectively defined by formulas (1.3) and (1.4). To prove this theorem, one relies on known formulas (see , ) for the inverse Fourier transforms of the functions $`|\xi |^{1/3}J_{1/3}(|\xi |),`$ $`|\xi |^{1/3}J_{1/3}(|\xi |),`$ and $`|\xi |^{1/3}K_{1/3}(|\xi |),`$ formulas that need be generalized to the case $`n>1.`$

## 3 Fourier transforms involving Bessel functions

As we have indicated in the previous section, we are going to need explicit formulas for the inverse Fourier transforms of the following functions: $`|\xi |^{\pm \nu }J_\nu (|\xi |),`$ $`|\xi |^\nu K_\nu (|\xi |),`$ and $`|\xi |^\nu N_\nu (|\xi |),`$ where $`\xi =(\xi _1,\mathrm{},\xi _n),`$ $`n1.`$ Since, when $`n=1,`$ explicit formulas for the Fourier transforms are known and can be found, for example, in  and , we concentrate on the case $`n>1.`$

1. The inverse Fourier transform of $`|\xi |^\nu K_\nu (|\xi |).`$

Assume that $`|\text{Re}\nu |<1/2`$ and define

(3.1)
$$^1[|\xi |^\nu K_\nu (|\xi |)](x)=\frac{1}{(2\pi )^n}__ne^{ix,\xi }|\xi |^\nu K_\nu (|\xi |)𝑑\xi .$$

The assumption on the real part of $`\nu `$ secures convergence of the integral at the origin. On the other hand, in view of the asymptotic behavior of $`K_\nu (|\xi |)`$ for large values of $`|\xi |,`$ the integral converges absolutely. By introducing spherical coordinates, rewrite the integral on the right-hand of (3.1) as

$$_0^{\mathrm{}}r^{n+\nu 1}K_\nu (r)\{_{S_{n1}}e^{irx,\omega }𝑑\omega \}𝑑r.$$

Since

(3.2)
$$_{S_{n1}}e^{irx,\omega }𝑑\omega =\frac{(2\pi )^{\frac{n}{2}}}{|rx|^{\frac{n}{2}1}}J_{\frac{n}{2}1}(r|x|),$$

it follows that

(3.3)
$$^1[|\xi |^\nu K_\nu (|\xi |)](x)=\frac{|x|^{1\frac{n}{2}}}{(2\pi )^{\frac{n}{2}}}_0^{\mathrm{}}r^{\frac{n}{2}+\nu }J_{\frac{n}{2}1}(r|x|)K_\nu (r)𝑑r.$$

We now quote the following result found in Watson’s treatise :

###### Lemma 3.1.

If $`\text{Re}(\mu +1)>|\text{Re}\nu |`$ and $`\text{Re}b>|\text{Im}a|`$ then

(3.4)
$$_0^{\mathrm{}}t^{\mu +\nu +1}J_\mu (at)K_\nu (bt)𝑑t=\frac{(2a)^\mu (2b)^\nu \mathrm{\Gamma }(\mu +\nu +1)}{(a^2+b^2)^{\mu +\nu +1}}.$$

From this lemma, we obtain the following result

###### Theorem 3.1.

If $`|\text{Re}\nu |<1/2`$ then

(3.5)
$$^1[|\xi |^\nu K_\nu (|\xi |)](x)=\frac{2^{\nu 1}\mathrm{\Gamma }(\frac{n}{2}+\nu )}{\pi ^{\frac{n}{2}}}\left(1+|x|^2\right)^{\frac{n}{2}\nu }.$$

Proof. The integral in (3.3) is the same as the integral in (3.4), where $`\mu =n/21,`$ $`a=|x|,`$ and $`b=1.`$ Clearly the conditions of the lemma are satisfied and Theorem 3.1 follows at once. $`\mathrm{}`$

We remark that, when $`n=1,`$ (3.5) becomes

$$^1[|\xi |^\nu K_\nu (|\xi |)](x)=\frac{2^{\nu 1}\mathrm{\Gamma }(\frac{1}{2}+\nu )}{\pi ^{\frac{1}{2}}}\left(1+|x|^2\right)^{\frac{1}{2}\nu },$$

a formula found in  and .

2. Inverse Fourier transforms of $`|\xi |^\nu J_\nu (|\xi |)`$ and $`|\xi |^\nu J_\nu (|\xi |).`$

As before, assume that $`|\text{Re}\nu |<1/2.`$ Formally, the inverse Fourier transform of $`|\xi |^\nu J_\nu (|\xi |)`$ is

$$^1[|\xi |^\nu J_\nu (|\xi |)](x)=\frac{1}{(2\pi )^n}__ne^{ix,\xi }|\xi |^\nu J_\nu (|\xi |)𝑑\xi .$$

In general, the integral diverges at $`\mathrm{}`$ and so we introduce a converging factor $`e^{ϵ|\xi |}`$ and take a limit as $`ϵ0.`$ Since $`|\xi |^\nu J_\nu (|\xi |)`$ is locally integrable, it defines a tempered distribution and the limit exists in $`𝒮^{}(_n),`$ the space of tempered distributions on $`_n.`$ Thus the precise meaning of the inverse Fourier transform is

(3.6)
$$^1[|\xi |^\nu J_\nu (|\xi |)](x)=\frac{1}{(2\pi )^n}\underset{ϵ0}{lim}__ne^{ix,\xi ϵ|\xi |}|\xi |^\nu J_\nu (|\xi |)𝑑\xi .$$

By introducing spherical coordinates, one can see that the integral on the right-hand can be written as

$$_0^{\mathrm{}}e^{ϵr}r^{n+\nu 1}J_\nu (r)\{_{S_{n1}}e^{irx,\omega }𝑑\omega \}𝑑r.$$

Since

(3.7)
$$_{S_{n1}}e^{irx,\omega }𝑑\omega =\frac{(2\pi )^{\frac{n}{2}}}{|rx|^{\frac{n}{2}1}}J_{\frac{n}{2}1}(r|x|),$$

it follows that

(3.8)
$$^1[|\xi |^\nu J_\nu (|\xi |)](x)=\frac{|x|^{1\frac{n}{2}}}{(2\pi )^{\frac{n}{2}}}\underset{ϵ0}{lim}_0^{\mathrm{}}e^{ϵr}r^{\frac{n}{2}+\nu }J_{\frac{n}{2}1}(r|x|)J_\nu (r)𝑑r.$$

In a similar manner we also have

(3.9)
$$^1[|\xi |^\nu J_\nu (|\xi |)](x)=\frac{|x|^{1\frac{n}{2}}}{(2\pi )^{\frac{n}{2}}}\underset{ϵ0}{lim}_0^{\mathrm{}}e^{ϵr}r^{\frac{n}{2}\nu }J_{\frac{n}{2}1}(r|x|)J_\nu (r)𝑑r.$$

Both integrals appearing in formulas (3.8) and (3.9) are particular cases of the following integral

(3.10)
$$I_ϵ(a,b)=_0^{\mathrm{}}e^{ϵt}t^\lambda J_\mu (at)J_\nu (bt)𝑑t,$$

where $`a`$ and $`b`$ are positive real numbers, $`\lambda ,`$ $`\mu ,`$ and $`\nu ,`$ complex numbers such that $`\text{Re}(\mu +\nu +1)>\text{Re}(\lambda ).`$ This integral is a variant of the discontinuous integral of Weber and Schafheitlin studied by Watson in his treatise on Bessel functions .

###### Lemma 3.2.

As $`ϵ0`$ $`I_ϵ(a,b)`$ tends, in the sense of distributions, to either one of the following limits:

(3.11)
$$\frac{b^\nu \mathrm{\Gamma }(\frac{\mu +\nu \lambda +1}{2})}{2^\lambda a^{\nu \lambda +1}\mathrm{\Gamma }(\nu +1)\mathrm{\Gamma }(\frac{\lambda +\mu \nu +1}{2})}F(\frac{\mu +\nu \lambda +1}{2},\frac{\nu \lambda \mu +1}{2};\nu +1;\frac{b^2}{a^2}),$$

if $`0<b<a,`$ or

(3.12)
$$\frac{a^\mu \mathrm{\Gamma }(\frac{\mu +\nu \lambda +1}{2})}{2^\lambda b^{\mu \lambda +1}\mathrm{\Gamma }(\mu +1)\mathrm{\Gamma }(\frac{\lambda +\nu \mu +1}{2})}F(\frac{\mu +\nu \lambda +1}{2},\frac{\mu \lambda \nu +1}{2};\mu +1;\frac{a^2}{b^2}),$$

if $`0<a<b.`$

Proof. We adapt to our situation the proof of the Weber-Schafheitlin theorem as found in section 13.4 of Watson’s treatise . It consists of expanding the integrand in (3.10) in power series of $`b,`$ and passing to the limit as $`ϵ0.`$

1. Consider the case when $`0<b<a.`$ If we replace $`b`$ by $`z,`$ the integral (3.10) is an analytic function of $`z`$ when $`\text{Re}z>0`$ and $`|\text{Im}z|<ϵ.`$ Introduce new constants $`\alpha ,`$ $`\beta ,`$ and $`\gamma `$ defined by

$$2\alpha =\mu +\nu \lambda +1,2\beta =\nu \lambda \mu +1,\gamma =\nu +1$$

or, equivalently,

$$\lambda =\gamma \alpha \beta ,\mu =\alpha \beta ,\nu =\gamma 1,$$

and rewrite (3.10) as follows:

$`I_ϵ(a,z)={\displaystyle _0^{\mathrm{}}}e^{ϵt}t^{\alpha +\beta \gamma }J_{\alpha \beta }(at)J_{\gamma 1}(zt)𝑑t`$
$`=`$ $`{\displaystyle _0^{\mathrm{}}}e^{ϵt}J_{\alpha \beta }(at)\left\{{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^m(z/2)^{\gamma +2m1}t^{\alpha +\beta +2m1}}{m!\mathrm{\Gamma }(\gamma +m)}}\right\}𝑑t`$
$`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^m(z/2)^{\gamma +2m1}}{m!\mathrm{\Gamma }(\gamma +m)}}{\displaystyle _0^{\mathrm{}}}e^{ϵt}J_{\alpha \beta }(at)t^{\alpha +\beta +2m1}𝑑t.`$
The interchange betwen the integration and summation signs is justified because, when $`|z|<ϵ,`$ the series

$$\underset{m=0}{\overset{\mathrm{}}{}}\frac{(1)^m(z/2)^{\gamma +2m1}}{m!\mathrm{\Gamma }(\gamma +m)}_0^{\mathrm{}}e^{ϵt}|J_{\alpha \beta }(at)t^{\alpha +\beta +2m1}|𝑑t$$

is absolutely convergent.

2. We now evaluate the last integral on the right-hand side of (3). By expanding $`J_{\alpha \beta }(at)`$ in power series and integrating term by term we obtain

$`{\displaystyle _0^{\mathrm{}}}e^{ϵt}J_{\alpha \beta }(at)t^{\alpha +\beta +2m1}𝑑t`$
$`={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^k(a/2)^{\alpha \beta +2k}}{k!\mathrm{\Gamma }(\alpha \beta +k+1)}}{\displaystyle _0^{\mathrm{}}}e^{ϵt}t^{2\alpha +2(m+k)1}𝑑t`$
$`={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^k(a/2)^{\alpha \beta +2k}}{k!\mathrm{\Gamma }(\alpha \beta +k+1)}}{\displaystyle \frac{\mathrm{\Gamma }(2\alpha +2m+2k)}{ϵ^{2\alpha +2m+2k}}}.`$
By using the duplication formula $`\mathrm{\Gamma }(2z)=2^{2z1}\pi ^{1/2}\mathrm{\Gamma }(z)\mathrm{\Gamma }(z+\frac{1}{2})`$ together with the notation $`(z)_k=z(z+1)\mathrm{}(z+k1)=\mathrm{\Gamma }(z+k)/\mathrm{\Gamma }(z)`$ we rewrite the last expression as:

$`{\displaystyle _0^{\mathrm{}}}e^{ϵt}J_{\alpha \beta }(at)t^{\alpha +\beta +2m1}𝑑t`$
$`=`$ $`{\displaystyle \frac{(a/2)^{\alpha \beta }}{(ϵ^2)^{\alpha +m}}}{\displaystyle \frac{\mathrm{\Gamma }(2\alpha +2m)}{\mathrm{\Gamma }(\alpha \beta +1)}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\alpha +m)_k(\alpha +m+\frac{1}{2})_k}{k!(\alpha \beta +1)_k}}\left({\displaystyle \frac{a^2}{c^2}}\right)^k`$
$`=`$ $`{\displaystyle \frac{(a/2)^{\alpha \beta }}{(ϵ^2)^{\alpha +m}}}{\displaystyle \frac{\mathrm{\Gamma }(2\alpha +2m)}{\mathrm{\Gamma }(\alpha \beta +1)}}F(\alpha +m,\alpha +m+{\displaystyle \frac{1}{2}};\alpha \beta +1;{\displaystyle \frac{a^2}{c^2}}).`$
By using formula (6.16) in the Appendix II:

$$F(a,b;c;z)=(1z)^aF(a,cb;c;\frac{z}{z1}),$$

we obtain that

$`{\displaystyle _0^{\mathrm{}}}e^{ϵt}J_{\alpha \beta }(at)t^{\alpha +\beta +2m1}𝑑t=`$
$`=`$ $`{\displaystyle \frac{(a/2)^{\alpha \beta }}{(a^2+ϵ^2)^{\alpha +m}}}{\displaystyle \frac{\mathrm{\Gamma }(2\alpha +2m)}{\mathrm{\Gamma }(\alpha \beta +1)}}F(\alpha +m,\alpha +m+{\displaystyle \frac{1}{2}};\alpha \beta +1;{\displaystyle \frac{a^2}{a^2+ϵ^2}}).`$
Substituting (3) into (3) we get

(3.15) $`I_ϵ(a,z)=`$
$$=\underset{m=0}{\overset{\mathrm{}}{}}\frac{(1)^m(z/2)^{\gamma +2m1}(a/2)^{\alpha \beta }\mathrm{\Gamma }(2\alpha +2m)}{m!\mathrm{\Gamma }(\gamma +m)(a^2+ϵ^2)^{\alpha +m}\mathrm{\Gamma }(\alpha \beta +1)}F(\alpha +m,\alpha +m+\frac{1}{2};\alpha \beta +1;\frac{a^2}{a^2+ϵ^2}),$$

whenever $`|z|<ϵ.`$

3. Following , one can show that (3.15) is valid provided that $`z`$ satisfies the conditions

$$\text{Re}(z)>0,|\text{Im}(z)|<ϵ,|z|<\sqrt{a^2+ϵ^2}ϵ.$$

Let $`\delta >0`$ be small enough so that $`b<\sqrt{a^2+\delta ^2}\delta ,`$ and take $`0<ϵ\delta `$ so that we also have $`b<\sqrt{a^2+ϵ^2}ϵ`$. In (3.15) we may now let $`z=b`$ and, when this is done, one can show, by the method of majorants, that the resulting series converges uniformly with respect to $`ϵ,`$ $`0<ϵ\delta ,`$ and therefore, as $`ϵ0,`$ the limit of the series is equal to its value at $`ϵ=0.`$ Thus

(3.16) $`\underset{ϵ0}{lim}I_ϵ(a,b)=`$
$$=\underset{m=0}{\overset{\mathrm{}}{}}\frac{(1)^m(b/2)^{\gamma +2m1}(a/2)^{\alpha \beta }\mathrm{\Gamma }(2\alpha +2m)}{m!\mathrm{\Gamma }(\gamma +m)a^{2\alpha +2m}\mathrm{\Gamma }(\alpha \beta +1)}F(\alpha +m,\alpha +m+\frac{1}{2};\alpha \beta +1;1).$$

By formula (6.15) in the Appendix II, we have

$$F(\alpha +m,\frac{1}{2}\beta m;\alpha \beta +1;1)=\frac{\mathrm{\Gamma }(\alpha \beta +1)\mathrm{\Gamma }(1/2)}{\mathrm{\Gamma }(1\beta m)\mathrm{\Gamma }(\alpha +m+\frac{1}{2})}.$$

On the other hand,

$$\mathrm{\Gamma }(1/2)\mathrm{\Gamma }(2\alpha +2m)=2^{2\alpha +2m1}\mathrm{\Gamma }(\alpha +m)\mathrm{\Gamma }(\alpha +m+\frac{1}{2}),$$

$$\mathrm{\Gamma }(1\beta m)=\frac{(1)^m\pi \mathrm{csc}(\pi \beta )}{\mathrm{\Gamma }(\beta +m)},$$

$$\mathrm{\Gamma }(1\beta )\mathrm{\Gamma }(\beta )=\pi \mathrm{csc}(\pi \beta ),$$

therefore

$$\frac{\mathrm{\Gamma }(2\alpha +2m)\mathrm{\Gamma }(1/2)}{\mathrm{\Gamma }(\gamma +m)\mathrm{\Gamma }(1\beta m)\mathrm{\Gamma }(\alpha +m+\frac{1}{2})}=\frac{2^{2\alpha +2m1}\mathrm{\Gamma }(\alpha +m)\mathrm{\Gamma }(\beta +m)}{(1)^m\mathrm{\Gamma }(1\beta )\mathrm{\Gamma }(\beta )\mathrm{\Gamma }(\gamma +m)}.$$

Substituting these formulas into (3.16), we get

$`\underset{ϵ0}{lim}I_ϵ(a,b)`$
$`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{b^{\gamma 1}\mathrm{\Gamma }(\alpha +m)\mathrm{\Gamma }(\beta +m)}{2^{\gamma \alpha \beta }a^{\alpha +\beta }m!\mathrm{\Gamma }(1\beta )\mathrm{\Gamma }(\beta )\mathrm{\Gamma }(\gamma +m)}}\left({\displaystyle \frac{b^2}{a^2}}\right)^m`$
$`=`$ $`{\displaystyle \frac{b^{\gamma 1}\mathrm{\Gamma }(\alpha )}{2^{\gamma \alpha \beta }a^{\alpha +\beta }\mathrm{\Gamma }(\gamma )\mathrm{\Gamma }(1\beta )}}F(\alpha ,\beta ;\gamma ;{\displaystyle \frac{b^2}{a^2}}).`$
Finally, returning to the constants $`\mu ,`$ $`\nu ,`$ and $`\lambda ,`$ we obtain the expression (3.11) in the first part of the lemma.

4. In the case $`0<a<b,`$ we proceed in an analogous manner to obtain the expression (3.12) in the second part of the lemma. $`\mathrm{}`$

We now use Lemma 3.2 to evaluate the inverse Fourier transforms of $`|\xi |^\nu J_\nu (|\xi |)`$ and $`|\xi |^\nu J_\nu (|\xi |)`$ respectively defined by formulas (3.8) and (3.9).

###### Theorem 3.2.

The inverse Fourier transform of $`|\xi |^\nu J_\nu (|\xi |)`$ is the distribution defined by

(3.22) $`^1[|\xi |^\nu J_\nu (|\xi |)](x)=`$
$`=\{\begin{array}{cc}\mathrm{sin}({\displaystyle \frac{n\pi }{2}}){\displaystyle \frac{2^\nu \mathrm{\Gamma }(\frac{n}{2}+\nu )}{\pi ^{\frac{n}{2}+1}}}(1|x|^2)^{\frac{n}{2}\nu },\hfill & 0<|x|<1,\hfill \\ & \\ \mathrm{sin}(\nu \pi ){\displaystyle \frac{2^\nu \mathrm{\Gamma }(\frac{n}{2}+\nu )}{\pi ^{\frac{n}{2}+1}}}(|x|^21)^{\frac{n}{2}\nu },\hfill & 1<|x|.\hfill \end{array}`$
Proof. The integral on the right-hand side of (3.8) corresponds to the integral in (3.10) where

$$\mu =\frac{n}{2}1,\nu =\nu ,\lambda =\frac{n}{2}\nu ,a=|x|,b=1.$$

If $`1<|x|,`$ then from formula (3.11) we obtain for the limit in (3.8):

$`\underset{ϵ0}{lim}{\displaystyle _0^{\mathrm{}}}e^{ϵr}r^{\frac{n}{2}+\nu }J_{\frac{n}{2}1}(r|x|)J_\nu (r)𝑑r=`$
$`={\displaystyle \frac{\mathrm{\Gamma }(\frac{n}{2}+\nu )}{2^{\frac{n}{2}\nu }|x|^{\frac{n}{2}+2\nu +1}\mathrm{\Gamma }(\nu +1)\mathrm{\Gamma }(\nu )}}F({\displaystyle \frac{n}{2}}+\nu ,\nu +1;\nu +1;{\displaystyle \frac{1}{|x|^2}}).`$
From the known relation for hypergeometric series

(3.23)
$$F(a,b;c;z)=(1z)^{cab}F(ca,cb;c;z)$$

we have that

$$F(\frac{n}{2}+\nu ,\nu +1;\nu +1;\frac{1}{|x|^2})=\left(\frac{|x|^21}{|x|^2}\right)^{\frac{n}{2}\nu }.$$

On the other hand

$$\mathrm{\Gamma }(\nu +1)\mathrm{\Gamma }(\nu )=\frac{\pi }{\mathrm{sin}(\pi \nu )}.$$

Thus, for $`1<|x|,`$ the above limit is equal to

$`\underset{ϵ0}{lim}{\displaystyle _0^{\mathrm{}}}e^{ϵr}r^{\frac{n}{2}+\nu }J_{\frac{n}{2}1}(r|x|)J_\nu (r)𝑑r=`$
$`=\mathrm{sin}(\pi \nu ){\displaystyle \frac{|x|^{\frac{n}{2}1}\mathrm{\Gamma }(\frac{n}{2}+\nu )}{2^{\frac{n}{2}\nu }\pi }}\left(|x|^21\right)^{\frac{n}{2}\nu },`$
If $`0<|x|<1,`$ then formula (3.12) yields

$`\underset{ϵ0}{lim}{\displaystyle _0^{\mathrm{}}}e^{ϵr}r^{\frac{n}{2}+\nu }J_{\frac{n}{2}1}(r|x|)J_\nu (r)𝑑r=`$
$`={\displaystyle \frac{|x|^{\frac{n}{2}1}\mathrm{\Gamma }(\frac{n}{2}+\nu )}{2^{\frac{n}{2}\nu }\mathrm{\Gamma }(\frac{n}{2})\mathrm{\Gamma }(1\frac{n}{2})}}F({\displaystyle \frac{n}{2}}+\nu ,{\displaystyle \frac{n}{2}};{\displaystyle \frac{n}{2}};|x|^2)`$
$`=\mathrm{sin}({\displaystyle \frac{n\pi }{2}}){\displaystyle \frac{|x|^{\frac{n}{2}1}\mathrm{\Gamma }(\frac{n}{2}+\nu )}{2^{\frac{n}{2}\nu }\pi }}\left(1|x|^2\right)^{\frac{n}{2}\nu },`$
after making the replacements

$$\mathrm{\Gamma }(\frac{n}{2})\mathrm{\Gamma }(1\frac{n}{2})=\frac{\pi }{\mathrm{sin}(n\pi /2)}$$

and

$$F(\frac{n}{2}+\nu ,\frac{n}{2};\frac{n}{2};|x|^2)=(1|x|^2)^{\frac{n}{2}\nu }.$$

The last expression results from (3.23). Substituting (3) and (3) into (3.8) we obtain (3.22) and the theorem is proved $`\mathrm{}`$

By reasoning in the same manner the following result, whose proof is left to the reader, holds.

###### Theorem 3.3.

The inverse Fourier transform of $`|\xi |^\nu J_\nu (|\xi |)`$ is the distribution defined by

$$^1[|\xi |^\nu J_\nu (|\xi |)](x)=\{\begin{array}{cc}\frac{(1|x|^2)^{\nu \frac{n}{2}}}{2^\nu \pi ^{\frac{n}{2}}\mathrm{\Gamma }(\nu \frac{n}{2}+1)},\hfill & 0<|x|<1,\hfill \\ & \\ 0,\hfill & 1<|x|.\hfill \end{array}$$

Remark. From theorem 3.3 we immediately derive a formula for the inverse Fourier transform of $`|\xi |^\nu J_\nu (|\xi |),`$ namely,

$$^1[|\xi |^\nu J_\nu (|\xi |)](x)=\{\begin{array}{cc}\frac{(1|x|^2)^{\nu \frac{n}{2}}}{2^\nu \pi ^{\frac{n}{2}}\mathrm{\Gamma }(1\nu \frac{n}{2})},\hfill & 0<|x|<1,\hfill \\ & \\ 0,\hfill & 1<|x|.\hfill \end{array}$$

Since

$$\mathrm{\Gamma }(1\nu \frac{n}{2})\mathrm{\Gamma }(\nu +\frac{n}{2})=\frac{\pi }{\mathrm{sin}[(\nu +\frac{n}{2})\pi ]},$$

we may rewrite (3) as

(3.26) $`^1[|\xi |^\nu J_\nu (|\xi |)](x)=`$
(3.30) $`=`$ $`\{\begin{array}{cc}{\displaystyle \frac{2^\nu \mathrm{sin}[(\nu +\frac{n}{2})\pi ]\mathrm{\Gamma }(\nu +\frac{n}{2})}{\pi ^{\frac{n}{2}+1}}}(1|x|^2)^{\nu \frac{n}{2}},\hfill & 0<|x|<1,\hfill \\ & \\ 0,\hfill & 1<|x|.\hfill \end{array}`$
3. The inverse Fourier transform of $`|\xi |^\nu N_\nu (|\xi |).`$

Theorems 3.2 and 3.3 imply the following result

###### Theorem 3.4.

The inverse Fourier transform of $`|\xi |^\nu N_\nu (|\xi |)`$ is the distribution defined by

(3.35) $`^1[|\xi |^\nu N_\nu (|\xi |)](x)=`$
$`=\{\begin{array}{cc}\mathrm{cos}({\displaystyle \frac{n\pi }{2}}){\displaystyle \frac{2^\nu \mathrm{\Gamma }(\frac{n}{2}+\nu )}{\pi ^{\frac{n}{2}+1}}}(1|x|^2)^{\frac{n}{2}\nu },\hfill & 0<|x|<1,\hfill \\ & \\ \mathrm{cos}(\nu \pi ){\displaystyle \frac{2^\nu \mathrm{\Gamma }(\frac{n}{2}+\nu )}{\pi ^{\frac{n}{2}+1}}}(|x|^21)^{\frac{n}{2}\nu },\hfill & 1<|x|.\hfill \end{array}`$
Proof. Recall that

$$N_\nu (|\xi |)=\frac{J_\nu (|\xi |)\mathrm{cos}(\nu \pi )J_\nu (|\xi |)}{\mathrm{sin}(\nu \pi )}.$$

Thus, from formulas (3.22) and (3.26) we get

$`^1[|\xi |^\nu N_\nu (|\xi |)](x)=`$
$`=\{\begin{array}{cc}{\displaystyle \frac{\mathrm{cos}(\nu \pi )}{\mathrm{sin}(\nu \pi )}}{\displaystyle \frac{2^\nu \mathrm{\Gamma }(\frac{n}{2}+\nu )}{\pi ^{\frac{n}{2}+1}}}(1|x|^2)^{\frac{n}{2}\nu },\hfill & 0<|x|<1,\hfill \\ & \\ \mathrm{cos}(\nu \pi ){\displaystyle \frac{2^\nu \mathrm{\Gamma }(\frac{n}{2}+\nu )}{\pi ^{\frac{n}{2}+1}}}(|x|^21)^{\frac{n}{2}\nu },\hfill & 1<|x|,\hfill \end{array}`$
$`\{\begin{array}{cc}{\displaystyle \frac{\mathrm{sin}[(\nu +\frac{n}{2})\pi ]}{\mathrm{sin}(\nu \pi )}}{\displaystyle \frac{2^\nu \mathrm{\Gamma }(\frac{n}{2}+\nu )}{\pi ^{\frac{n}{2}+1}}}(1|x|^2)^{\frac{n}{2}\nu },\hfill & 0<|x|<1,\hfill \\ & \\ 0,\hfill & 1<|x|.\hfill \end{array}`$
Now

$$\mathrm{sin}[(\nu +\frac{\pi }{2})]\pi =\mathrm{sin}(\nu \pi )\mathrm{cos}(\frac{n\pi }{2})+\mathrm{cos}(\nu \pi )\mathrm{sin}(\frac{n\pi }{2})$$

and so formula (3.35) follows at once. $`\mathrm{}`$

## 4 A fundamental solution with support in $`\overline{D_{}^n}`$

We return to the problem (2.1) and, in order to obtain a formula for a fundamental solution with support in the hyperbolic region, we have to modify formula (2.13). Consider the function

(4.1)
$$\stackrel{~}{F}_{}(\xi ,y)=\{\begin{array}{cc}\frac{3^{2/3}\mathrm{\Gamma }(4/3)}{2^{1/3}}(\frac{t}{|\xi |})^{1/3}J_{1/3}(t|\xi |)\hfill & \text{for }y0\hfill \\ & \\ 0\hfill & \text{for }y0,\hfill \end{array}$$

where $`t=(2/3)(y)^{3/2}.`$ From formula (6.10) in the Appendix I, it follows that the limit of $`\stackrel{~}{F}_{}(\xi ,y)`$ as $`y0`$ is equal to zero and so condition (i) in Section 2 holds. Also it follows from the same formula that

$$_y\{(\frac{t}{|\xi |})^{1/3}J_{1/3}(t|\xi |)\}_{|y=0}=\frac{2^{1/3}}{3^{2/3}\mathrm{\Gamma }(4/3)},$$

hence, the $`y`$-derivative of $`\stackrel{~}{F}_{}(\xi ,y)`$ at $`y=0`$ is equal to $`1,`$ and so condition (ii) is also satisfied. By calculating the inverse Fourier transform of $`\stackrel{~}{F}_{}(\xi ,y)`$ we have the following result:

###### Theorem 4.1.

The inverse Fourier transform of $`\stackrel{~}{F}_{}(\xi ,y)`$ is the distribution

(4.2)
$$F_{}(x,y)=\{\begin{array}{cc}\frac{3^n\mathrm{\Gamma }(\frac{4}{3})}{2^{2/3}\pi ^{n/2}\mathrm{\Gamma }(\frac{4}{3}\frac{n}{2})}|9|x|^2+4y^3|^{\frac{1}{3}\frac{n}{2}}\hfill & \text{in }D_{}^n\hfill \\ & \\ 0\hfill & \text{elsewhere,}\hfill \end{array}$$

where $`D_{}^n=\{(x,y)^{n+1}:9|x|^2+4y^3<0\}.`$ It is a fundamental solution for the operator (1.1) whose support is the closure of the region $`D_{}^n.`$

Proof. We must evaluate the inverse Fourier transform of $`(t/|\xi |)^{1/3}J_{1/3}(t|\xi |).`$ Recall that if $`G(x)`$ is the inverse Fourier transform of $`f(\xi ),`$ then $`(1/a^n)G(x/a)`$ is the inverse Fourier transform of $`f(a\xi ).`$ By appying formula (3.3) for $`\nu =1/3,`$ we obtain

$`^1[(t/|\xi |)^{1/3}J_{1/3}(t|\xi |)](x)=`$
$`=\{\begin{array}{cc}{\displaystyle \frac{t^{\frac{2}{3}n}}{2^{1/3}\pi ^{\frac{n}{2}}\mathrm{\Gamma }(\frac{4}{3}\frac{n}{2})}}\left(1{\displaystyle \frac{|x|^2}{t^2}}\right)^{\frac{1}{3}\frac{n}{2}},\hfill & 0<|x|<t,\hfill \\ & \\ 0,\hfill & t<|x|.\hfill \end{array}`$
$`=\{\begin{array}{cc}{\displaystyle \frac{1}{2^{1/3}\pi ^{\frac{n}{2}}\mathrm{\Gamma }(\frac{4}{3}\frac{n}{2})}}\left(t^2|x|^2\right)^{\frac{1}{3}\frac{n}{2}},\hfill & 0<|x|<t,\hfill \\ & \\ 0,\hfill & t<|x|.\hfill \end{array}`$
$`=\{\begin{array}{cc}{\displaystyle \frac{3^{n\frac{2}{3}}}{2^{1/3}\pi ^{\frac{n}{2}}\mathrm{\Gamma }(\frac{4}{3}\frac{n}{2})}}\left|9|x|^2+4y^3\right|^{\frac{1}{3}\frac{n}{2}}\hfill & \text{in}D_{}^n\hfill \\ & \\ 0,\hfill & \text{elsewhere}.\hfill \end{array}`$
Multiplication by the constant $`3^{2/3}\mathrm{\Gamma }(4/3)/2^{1/3}`$ in formula (4.1) yields (4.2) which proves the theorem. $`\mathrm{}`$

Remarks 1. If $`n=1,`$ then $`F_{}(x,y)`$ coincides with the distribution defined by formula (1.4).

2. For all values of $`n`$ the support of $`F_{}(x,y)`$ is the closure of the region $`D_{}^n.`$ This follows from results in  about the generalized function $`P_+^\lambda `$ where $`P(x)`$ is the quadratic polynomial

$$P(x)=x_1^2+\mathrm{}+x_p^2x_{p+1}^2\mathrm{}x_{p+q}^2,$$

$`p,q1.`$ See also the example in  about the distribution $`\text{Pf}.s^\lambda ,`$ where $`s`$ denotes the hyperbolic distance in $`^N,`$ $`N>1.`$

## 5 A fundamental solution with support in $`\overline{D_+^n}`$

Consider the function $`\stackrel{~}{F}_+(\xi ,y)`$ defined by

(5.1)
$$\stackrel{~}{F}_+(\xi ,y)=\{\begin{array}{cc}\gamma (\frac{s}{|\xi |})^{1/3}K_{1/3}(s|\xi |)\hfill & \text{if }y0\hfill \\ \delta (\frac{t}{|\xi |})^{1/3}N_{1/3}(t|\xi |)\hfill & \text{if }y0,\hfill \end{array}$$

where $`s=(2/3)y^{3/2},`$ $`t=(2/3)(y)^{3/2},`$ and the constants $`\gamma `$ and $`\delta `$ are respectively given by

(5.2)
$$\gamma =\frac{2^{2/3}}{3^{4/3}\mathrm{\Gamma }(2/3)}\text{and}\delta =\frac{2\pi }{2^{1/3}3^{4/3}\mathrm{\Gamma }(2/3)},$$

and where $`N_{1/3}`$ is the Neumann function defined by (6.4). By using the formulas in the Appendix I, it is a matter of verification that the conditions (i) and (ii) in Section 2 are satisfied. Thus its inverse Fourier transform, denoted by $`F^{\mathrm{}}(x,y),`$ defines a fundamental solution for the operator (1.1). The following theorem gives an explicit expression for $`F^{\mathrm{}}(x,y).`$

###### Theorem 5.1.

The inverse Fourier transform of $`\stackrel{~}{F}_+(\xi ,y)`$ is the distribution defined by

(5.3)
$$F^{\mathrm{}}(x,y)=\{\begin{array}{cc}\frac{3^{n2}}{2^{2/3}\pi ^{n/2}}\frac{\mathrm{\Gamma }(\frac{n}{2}\frac{1}{3})}{\mathrm{\Gamma }(2/3)}(9|x|^2+4y^3)^{\frac{1}{3}\frac{n}{2}}\hfill & \text{in }D_+^n\hfill \\ & \\ \mathrm{cos}(\frac{n\pi }{2})\frac{2^{1/3}3^{n2}}{\pi ^{n/2}}\frac{\mathrm{\Gamma }(\frac{n}{2}\frac{1}{3})}{\mathrm{\Gamma }(2/3)}|9|x|^2+4y^3|^{\frac{1}{3}\frac{n}{2}}\hfill & \text{in }D_{}^n\text{.}\hfill \end{array}$$

Proof. 1. We start by evaluating the inverse Fourier transform of $`(s/|\xi |)^{1/3}K_{1/3}(s|\xi |),`$ $`s=2y^{3/2}/3,`$ $`y0.`$ Since $`K_\nu (z)=K_\nu (z),`$ it follows from formula (3.5) that

$$^1[(s/|\xi |)^{1/3}K_{1/3}(s|\xi |)]=\frac{2^{\frac{4}{3}}\mathrm{\Gamma }(\frac{n}{2}\frac{1}{3})}{\pi ^{n/2}}\left(s^2+|x|^2\right)^{\frac{1}{3}\frac{n}{2}}.$$

After multiplying by the constant $`\gamma `$ in (5.2) and reintroducing the variables $`x`$ and $`y,`$ we can see that the right-hand side of the last expression is equal to

(5.4)
$$\frac{3^{n2}\mathrm{\Gamma }(\frac{n}{2}\frac{1}{3})}{2^{2/3}\pi ^{n/2}\mathrm{\Gamma }(2/3)}\left(9|x|^2+4y^3\right)^{\frac{1}{3}\frac{n}{2}},y0,$$

that is $`F^{\mathrm{}}(x,y)`$ in $`D_+^n\{y0\}.`$

2. Next we evaluate the inverse Fourier transform of $`(t/|\xi |)^{1/3}N_{1/3}(t|\xi |),`$ $`t=2(y)^{3/2}/3,`$ $`y0.`$ From formula (3.35) we obtain

$`^1[(t|\xi |)^{1/3}N_{1/3}(t|\xi |)]=`$
$`=\{\begin{array}{cc}\mathrm{cos}({\displaystyle \frac{n\pi }{2}}){\displaystyle \frac{2^{1/3}\mathrm{\Gamma }(\frac{n}{2}\frac{1}{3})}{\pi ^{\frac{n}{2}+1}}}(t^2|x|^2)^{\frac{1}{3}\frac{n}{2}},\hfill & 0<|x|<t,\hfill \\ & \\ \left({\displaystyle \frac{1}{2}}\right){\displaystyle \frac{2^{1/3}\mathrm{\Gamma }(\frac{n}{2}\frac{1}{3})}{\pi ^{\frac{n}{2}+1}}}(|x|^2t^2)^{\frac{1}{3}\frac{n}{2}},\hfill & t<|x|.\hfill \end{array}`$
After multiplying both sides by the constant $`\delta `$ in (5.2) and reverting to the variables $`x`$ and $`y,`$ we can see that the right-hand side of the last expression can be written as

$$\{\begin{array}{cc}\frac{3^{n2}\mathrm{\Gamma }(\frac{n}{2}\frac{1}{3})}{2^{2/3}\pi ^{n/2}\mathrm{\Gamma }(2/3)}\left(9|x|^2+4y^3\right)^{\frac{1}{3}\frac{n}{2}},\hfill & \text{in}D_+^n\{y0\}\hfill \\ & \\ \mathrm{cos}(\frac{n\pi }{2})\frac{2^{1/3}3^{n2}\mathrm{\Gamma }(\frac{n}{2}\frac{1}{3})}{\pi ^{n/2}\mathrm{\Gamma }(2/3)}\left|9|x|^2+4y^3\right|^{\frac{1}{3}\frac{n}{2}},\hfill & \text{in}D_{}^n,\hfill \end{array}$$

and so (5) coincides with $`F^{\mathrm{}}(x,y)`$ in $`D_+^n\{y0\}D_{}^n.`$ Therefore, from (5.4) and (5) we obtain (5.3) and the theorem is proved. $`\mathrm{}`$

Remarks 1. In view of the results in  relative to the generalized function $`P_+^\lambda `$ (see Remark 2 after the proof of Theorem 4.1) the supports of the distributions

$$(9|x|^2+4y^3)^{\frac{1}{3}\frac{n}{2}}\text{and}|9|x|^2+4y^3|^{\frac{1}{3}\frac{n}{2}}$$

that appear in formula (5.3) are the closures of $`D_+^n`$ and $`D_{}^n,`$ respectively. In other words, for no value of $`n`$ can these supports be just the boundaries of these regions.

2. If $`n`$ is odd, then $`\mathrm{cos}(n\pi /2)=0,`$ and we rewrite formula (5.3), using the notation $`F_+(x,y)`$ instead of $`F^{\mathrm{}}(x,y),`$ as follows:

(5.6)
$$F_+(x,y)=\{\begin{array}{cc}\frac{3^{n2}}{2^{2/3}\pi ^{n/2}}\frac{\mathrm{\Gamma }(\frac{n}{2}\frac{1}{3})}{\mathrm{\Gamma }(2/3)}(9|x|^2+4y^3)^{\frac{1}{3}\frac{n}{2}}\hfill & \text{in }D_+^n\hfill \\ & \\ 0\hfill & \text{elsewhere.}\hfill \end{array}$$

This is a fundamental solution whose support is $`\overline{D_+^n}.`$ In particular, if $`n=1`$ we obtain formula (1.3).

3. If $`n`$ is even, then $`F^{\mathrm{}}(x,y)`$ is not necessarily identically zero in $`D_{}^n`$ and its support may be the whole of $`^{n+1}.`$ Suppose that $`n=2k,k>0,`$ and let us compare the constant in formula (5.3), relative to the region $`D_{}^n,`$ to the constant in the expression (4.2) of $`F_{}(x,y).`$ Let

$$A=\frac{(1)^{k+1}2^{1/3}3^{2(k1)}}{\pi ^k}\frac{\mathrm{\Gamma }(k\frac{1}{3})}{\mathrm{\Gamma }(2/3)}$$

be the constant in (5.3) and let

$$B=\frac{3^{2k}}{2^{2/3}\pi ^k}\frac{\mathrm{\Gamma }(4/3)}{\mathrm{\Gamma }(\frac{4}{3}k)}=\frac{(1)^{k+1}3^{2k1}}{2^{2/3}\pi ^k}\frac{\mathrm{\Gamma }(k\frac{1}{3})}{\mathrm{\Gamma }(2/3)}$$

be the constant in (4.2). Since $`3A2B=0`$ it follows that the distribution

(5.7)
$$F_+(x,y)=3F^{\mathrm{}}(x,y)2F_{}(x,y)$$

is now a fundamental solution for the operator (1.1) supported in $`\overline{D_+^n}`$ and we have for this $`F_+(x,y)`$ the same expression as that of (5.6). In conclusion, for all values on $`n`$ we always get two fundamental solutions: one whose support is the closure of $`D_+^n`$ and another whose support is the closure of $`D_{}^n.`$

## 6 Appendix

I. Bessel functions

The function $`J_\nu (z)`$ of a complex variable $`z`$ defined by

(6.1)
$$J_\nu (z)=\underset{r=0}{\overset{\mathrm{}}{}}\frac{(1)^rz^{\nu +2r}}{2^{\nu +2r}r!\mathrm{\Gamma }(\nu +r+1)},|z|<\mathrm{},|\mathrm{arg}z|<\pi ,$$

is called the Bessel function of the first kind of order $`\nu `$.

We also need the Bessel functions $`I_\nu (z)`$ and $`K_\nu (z)`$ defined by

(6.2)
$$I_\nu (z)=\underset{r=0}{\overset{\mathrm{}}{}}\frac{z^{\nu +2r}}{2^{\nu +2r}r!\mathrm{\Gamma }(\nu +r+1)},|z|<\mathrm{},|\mathrm{arg}z|<\pi ,$$

and

(6.3)
$$K_\nu (z)=\frac{\pi \mathrm{csc}(\nu \pi )}{2}\{I_\nu (z)I_\nu (z)\},\nu 0,\pm 1,\pm 2,\mathrm{}$$

as well as the Neumann function

(6.4)
$$N_\nu (z)=\frac{J_\nu (z)\mathrm{cos}(\nu \pi )J_\nu (z)}{\mathrm{sin}(\nu \pi )}.$$

Note that throughout this work, we only deal with Bessel functions of order $`\pm 1/3`$. Recall that $`Ai(z)`$ was defined in formula (2.6) by

$$Ai(z)=\frac{z^{1/2}}{3}[I_{1/3}(\frac{2}{3}z^{3/2})I_{1/3}(\frac{2}{3}z^{3/2})]=\frac{1}{\pi }(\frac{z}{3})^{1/2}K_{1/3}(\frac{2}{3}z^{3/2}).$$

If we set $`s=\frac{2}{3}z^{3/2}`$, then we may rewrite $`Ai(z)`$ as

(6.5)
$$Ai(z)=\frac{1}{3^{2/3}2^{1/3}}s^{1/3}[I_{1/3}(s)I_{1/3}(s)].$$

From the series expansion of $`I_\nu (z)`$ it follows that

(6.6)
$$s^{1/3}I_{1/3}(s)=\frac{1}{2^{1/3}\mathrm{\Gamma }(2/3)}+\frac{s^2}{2^{5/3}\mathrm{\Gamma }(5/3)}+\mathrm{}$$

and

(6.7)
$$s^{1/3}I_{1/3}(s)=\frac{s^{2/3}}{2^{1/3}\mathrm{\Gamma }(4/3)}+\frac{s^{8/3}}{2^{7/3}\mathrm{\Gamma }(7/3)}+\mathrm{}$$

Consequently from (6.5), (6.6), and (6.7) we obtain

$$Ai(0)=\frac{3^{2/3}}{\mathrm{\Gamma }(2/3)},$$

the first expression in formula (2.8). Similarly, by differentiating $`Ai(z)`$ and setting $`z=0`$, we get the second expression in (2.8)

$$Ai^{}(0)=\frac{3^{4/3}}{\mathrm{\Gamma }(4/3)}.$$

In an analogous way, recall that $`Bi(z)`$ in (2.6) was defined by

$$Bi(z)=(\frac{z}{3})^{1/2}[I_{1/3}(\frac{2}{3}z^{3/2})+I_{1/3}(\frac{2}{3}z^{3/2})].$$

If we set $`s=\frac{2}{3}z^{3/2}`$, then the last expression becomes

(6.8)
$$Bi(z)=\frac{1}{2^{1/3}3^{1/6}}s^{1/3}[I_{1/3}(s)+I_{1/3}(s)],$$

and again from (6.6) and (6.7) we obtain

$$Bi(0)=\frac{3^{1/6}}{\mathrm{\Gamma }(2/3)},Bi^{}(0)=\frac{3^{5/6}}{\mathrm{\Gamma }(4/3)},$$

which are the two expressions in (2.9).

Finally, from these results we obtain the value of the Wronskian of $`Ai(z)`$ and $`Bi(z)`$ at $`z=0`$:

$$W(Ai(z),Bi(z))_{|_{z=0}}=\frac{23^{3/2}}{\mathrm{\Gamma }(2/3)\mathrm{\Gamma }(4/3)}=1/\pi ,$$

because $`\mathrm{\Gamma }(2/3)\mathrm{\Gamma }(4/3)=2\pi 3^{3/2}.`$

For future reference, we also need the following expansions derived from the series that defines $`J_\nu (z)`$:

(6.9)
$$t^{1/3}J_{1/3}(t)=\frac{1}{2^{1/3}\mathrm{\Gamma }(2/3)}\frac{t^2}{2^{5/3}\mathrm{\Gamma }(5/3)}+\mathrm{}$$

and

(6.10)
$$t^{1/3}J_{1/3}(t)=\frac{t^{2/3}}{2^{1/3}\mathrm{\Gamma }(4/3)}\frac{t^{8/3}}{2^{7/3}\mathrm{\Gamma }(7/3)}+\mathrm{}$$

II. Hypergeometric series and functions

Let $`a,`$ $`b,`$ and $`c`$ be arbitrary complex numbers and let $`z`$ be a complex variable. The power series

(6.11)
$$F(a,b;c;z)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{(a,n)(b,n)}{(c,n)}\frac{z^n}{n!},$$

where

(6.12)
$$(a,0)=1,(a,n)=\frac{\mathrm{\Gamma }(a+n)}{\mathrm{\Gamma }(a)}=a(a+1)\mathrm{}(a+n1),n=1,2,\mathrm{}$$

and we assume that $`c0,1,2,\mathrm{},`$ is called a hypergeometric series. It is known  that the series (6.11) is a solution, valid near $`z=0,`$ of the hypergeometric equation

(6.13)
$$z(1z)\frac{d^2u}{dz^2}+\{c(a+b+1)z\}\frac{du}{dz}abu=0,$$

for which every point is an ordinary point, except $`0,`$ $`1,`$ and $`\mathrm{},`$ that are regular singular points.

If either $`a`$ or $`b`$ is a negative integer, then the series (6.11) terminates; if $`c`$ is a negative integer, the series is meaningless because all terms after the $`(1c)`$ th have a zero denominator. As it is known , it is possible to redefine the series (6.11) so that it still is a solution of the hypergeometric equation. We exclude this possibility from our considerations because, in the cases that interest us, $`c`$ is never a negative integer.

The hypergeometric series is absolutely convergent for $`|z|<1`$ and so defines, in the open disk, an analytic function of $`z`$ which is regular at $`z=0.`$ The point $`z=1`$ is however a branch point and if a cut is made from $`1`$ to $`+\mathrm{}`$ along the $`x`$-axis, it can be shown  that series can be continued analytically and defines an analytic function throughout the cut plane that we still denote by $`F(x,b;c;z).`$ If $`\text{Re}(c)>\text{Re}(b)>0,`$ this analytic extension can be represented by Euler’s formula

(6.14)
$$F(a,b;c;z)=\frac{\mathrm{\Gamma }(c)}{\mathrm{\Gamma }(b)\mathrm{\Gamma }(cb)}_0^1t^{b1}(1b)^{cb1}(1tz)^a𝑑t,$$

for $`|\text{ }arg(1z)|<\pi .`$

In general, the hypergeometric series (6.11) diverges for $`|z|=1.`$ However, if $`\text{Re}(cab)>0,`$ we have absolute convergence for $`|z|=1.`$ Moreover,

(6.15)
$$F(a,b;c;1)=\frac{\mathrm{\Gamma }(c)\mathrm{\Gamma }(cab)}{\mathrm{\Gamma }(ca)\mathrm{\Gamma }(cb)}.$$

Hypergeometric functions satisfy among themselves quite a number of important relations of which we just list the following two that we used in Section 3:

(6.16)
$$F(a,b;c;z)=(1z)^aF(a,cb;c;\frac{z}{z1})$$

and

(6.17)
$$F(a,b;c;z)=(1z)^{cab}F(ca,cb;c;z).$$

For a complete list of such relations, the reader should consult Erdély .


text:
# Hydration of polyelectrolytes studied by molecular dynamics simulation

## 1 Introduction

Polyelectrolytes play an important role in industrial chemistry. The fields of application range from tailormade thickeners to paper finishing or ore preparation. Polysaccharide derivatives represent one interesting class of polyelectrolytes. In particular, cellulose products are important compounds. For our simulation study carboxy methyl cellulose (CMC) is chosen as an example for a polyelectrolyte derived from a natural polymer. Aqueous CMC solutions exhibit valuable properties, like a wide range of viscosity, non-toxicity and biodegradability. Particularly for the high-purity consumer-product market (cosmetics, food stuffs), CMC is used. However, pricing becomes more important in bulk applications (clay and ore treatment, oil-drilling). Hence it is desirable to replace some of the high-cost–high-selective chemicals with low-cost equivalents, like poly(acrylic acid) (PAA), which is the prototype of industrial synthetic polyelectrolytes. PAA is the other polymer studied in this work.

Most published work on aqueous CMC and PAA solutions was done experimentally using chromatography , 13C nuclear magnetic resonance  and rheological techniques . Theoretical approaches are scarce. We are aware of only one paper , which treats CMC by the worm-like-chain theory. This electrostatic theory successfully rationalizes some of the global properties of CMC, but as a rather generic approach it does not allow for detailed predictions on an atomistic time and length scale. Similar restrictions apply also to Monte Carlo simulations of polyelectrolyte chains in a cell model, where the solvent is treated as a dielectric continuum . Especially local interactions hydrogen bonds (hydrogen bridges) are neglected in theories and non-atomistic simulations. With two or three hydrogen-bond donor groups per repeat unit and even more acceptors sites (including charged COO- groups), this type of interaction is likely to be very important for the behaviour of CMC in water. Experimental techniques, on the other hand, suffer from different problems: NMR provides averaged local properties. Rheology derives and verifies scaling laws, but different polyelectrolytes lose their chemical “identity” and show the generic behaviour of excluded-volume chains. Atomistic molecular dynamics (MD) simulations cannot overcome all these problems, but they can provide some more detailed information.

The aim of this study is to understand better the structural and dynamic aspects of the hydration of CMC and PAA and to compare the two polymers. To this end, we investiagate both the chain properties and the interaction of chains with their immediate solvent environment. Atomistic simulation is confined to the study of small system sizes. However, in combination with coarse-graining methods  even some mesoscopic properties may be explored. Thus, a second goal of this study was to produce atomistic structural information, from which coarse grained models of e. g. PAA can be generated.

## 2 Computational Details

The exact molecular composition of CMC is a consequence of the molecule’s production history. Usually OH-groups are randomly substituted by carboxy methyl units  and the degree of substitution (DS) for most commercial CMC’s lies in the range DS=$`0.6`$$`1.2`$ substitutions per anhydroglucose unit (AGU). There is evidence for a preferential substitution during heterogeneous, wet conversion of cellulose to CMC: Heinze and coworkers used liquid chromatography to analyse the substitution pattern of Sodium-CMC. Their analysis was preceded by hydrolysis of the glucosidic link, by which information on the ring-substitution sequence is lost. For the AGU, the following substitution statistics was found : 23% of all O2-oxygens were substituted, followed by the O6 (15%) and O3 with 11% of all respective atoms. The glucose nomenclature is shown in figure 1. The diversity of molecular substitution patterns of CMC makes it difficult to model such a compound by an oligomer, because there is no “typical” molecule or structure.

To have at least some variety, two CMC oligomers were simulated. Substitution patterns were generated at random, with weights for the three exocyclic oxygens ($`P(O2)=0.23`$, $`P(O3)=0.12`$, $`P(O6)=0.15`$). They are summarized in table 1. The first molecule (CMC I) consists of seven AGUs and is substituted five times, the second one (CMC II) consist of eight AGUs with six OH groups substituted. The degrees of substitution are DS=$`0.70`$ and DS=$`0.75`$, respectively. Due to three sites, at which substitution can take place (neglecting the chain ends) there are eight different possible substitution patterns on a single AGU. Four of these (O2 O3, O3, O2, O6) are present in the first simulated CMC oligomer and one further pattern (O3 O6) in the second CMC oligomer.

All carboxylic acid groups were assumed to be dissociated, neutrality of the systems was maintained by five, respectively six sodium counterions. As the carboxylic acid group has a $`\mathrm{pK}_\mathrm{A}`$ of about $`4`$, the system is basic. Molecular weights are $`1440`$ amu (seven AGUs, CMC I) and $`1660`$ amu (eight AGUs, CMC II), respectively. Molecule CMC I was solvated in 3789, molecule II in 3917 water molecules. The concentrations are about $`2.1`$ and $`2.4`$ weight percent (excluding the sodium-counterions). Therefore both simulations were performed in the dilute regime (no added salt), where polyelectrolytes are expected to be stretched .

Polyacrylic acid (PAA) was modeled as one atactic oligomer strand, with $`23`$ repeat units ($`1636`$ amu) solvated in $`3684`$ water molecules (concentration excluding counterions $`2.4`$ weight percent). All COO- groups are negatively charged and there are $`23`$ sodium counterions. Runs were performed under the same conditions as for CMC, including the GROMOS force field  as detailed in table 2 and 3.

To ensure a consistent forcefield, we used simple point charge (SPC) water  with the GROMOS parameter set. This combination has proven useful in several simulations of sugars in water . Lennard-Jones interactions were truncated and shifted at a cutoff of $`r_c=0.9\mathrm{nm}`$. An isotropic pressure correction term was applied for site-site distances greater than the cutoff $`r_c`$. Electrostatics were dealt with using the Coulomb-Potential with a reaction-field correction :

$$V_q(r_{ij})=\frac{q_iq_j}{4\pi ϵϵ_0}\left(\frac{1}{r_{ij}}+\frac{ϵ_{\mathrm{RF}}1}{2ϵ_{\mathrm{RF}}+1}\frac{r_{ij}^2}{r_c^3}\right)$$

with the reaction field dielectric $`ϵ_{\mathrm{RF}}`$ equal to the value for water ($`ϵ_{\mathrm{RF}}=78.5`$). All atoms – except aliphatic hydrogens – are explicitly modeled. In the GROMOS forcefield, aliphatic hydrogens are accounted for by a change of the parameters for the parent carbons.

Bond lengths were held constant using the Shake procedure . Equations of motion were integrated using the leap-frog algorithm with a timestep of $`2`$ fs. Other simulation parameters were a neighbor list, which was updated every $`15`$ steps with a list cutoff of $`1\mathrm{nm}`$, weak coupling to a temperature ($`T=333.15\mathrm{K}`$) and pressure ($`p=1\mathrm{atm}`$) bath  with coupling times of $`\tau _T=0.5\mathrm{ps}`$ and $`\tau _p=2.5\mathrm{ps}`$ (water compressibility $`4.510^{10}\mathrm{kPa}^1`$). For later analysis, coordinates were written every 195 steps. The overall simulation time was 4.5 ns for CMC I and $`2.5`$ ns for CMC II, each after $`1`$ ns of equilibration.

We used different starting geometries for both CMC I and CMC II. The first one was started in a stretched conformation, the second, CMC II, was prepared in a folded starting conformation: During the first $`150\mathrm{ps}`$ of equilibration, a globule-like conformation was enforced by five artificial harmonic bonds (force constant of $`100\mathrm{kJ}\mathrm{mol}^1\mathrm{nm}^2`$ with a minimum energy at $`1\mathrm{nm}`$ between C1 atoms of different glucose rings). PAA was started in the stretched (all-trans) conformation, equilibrated for $`1`$ ns and simulated for about $`4.5`$ ns.

For analysis purpose, the center of mass of a repeat unit is caluclated using all its atoms, including any side groups if present. The CMC repeat unit begins at C1 and ends at O4 according to the numbering in figure 1. The acrylate repeat unit coincides with the propylate unit $`\text{-}\mathrm{CH}_2\text{-}\mathrm{CH}\text{-}\mathrm{COO}^{}`$.

The simulations were done with our molecular dynamics package YASP, which is described in ref. .

## 3 Results

### 3.1 Global Chain Properties

Carboxy methyl cellulose I (CMC I) and CMC II have very different structures: CMC I is an extended, almost straight chain, whereas CMC II assumes a cyclic conformation (Fig. 2). Poly(acrylic acid) (PAA) behaves more like CMC I, having mostly a straight conformation with only some bending (Fig. 3).

These structural features cause the end-to-end distance distributions (figure 4) of CMC I and CMC II to be very different. The end-to-end distance $`R_{\mathrm{ete}}`$ is calculated between the center of mass of the first and last AGU in the chain. For PAA it is calculated from the center of mass of the terminate repeat units. For CMC I there are two conformers, a major one at $`2.6\mathrm{nm}`$ and a second state at $`1.8\mathrm{nm}`$. This second state corresponds to a bent conformation, which shows up several times in the simulation. Its probability is about 15% of the stretched structure and on average this state is kept for about $`80\mathrm{ps}`$ until the oligomer’s end-to-end distance lengthens again. For CMC II the situation is different: The end-to-end distribution is very narrow. The short distance corresponds to a stable ring conformer. Poly(acrylic acid) shows a wider distribution than either CMC, resulting from the larger number of torsional degrees of freedom. The shoulder at $`2.8\mathrm{nm}`$ is due to the bent conformation shown in Fig. 3.

The shapes of the molecules are also reflected in the the eigenvalues of the gyration tensor (table 4). The stretched conformations of CMC I and PAA, lead to one large and two small eigenvalues, whereas CMC II has more isotropic eigenvalues, with only a factor of two between the largest and second largest value.

The correlation between the eigenvalues of the gyration tensor follows opposite trends (table 4): Whereas CMC I and PAA exhibit a negative correlation ($`0.86`$ and $`0.76`$) between the two largest eigenvalues, there is a small positive correlation ($`0.14`$) for CMC II. This shows how structural fluctuations happen: CMC I and PAA behave like bending rods – they shrink in one dimension and necessarily grow in another. In contrast the ring pulsation of CMC II, leads to the two largest eigenvalues growing and shrinking simultaneously.

Both the gyration radius and the end-to-end distance have multiple time regimes even in the first $`80\mathrm{ps}`$. At higher correlation times the correlation function’s statistics get worse.

The time correlation function of the radius of gyration (figure 5) of CMC I displays two time regimes: on the first $`15\mathrm{ps}`$ a fast one with a time constant of $`41\mathrm{ps}`$ followed by a slow decay with $`\tau =92\mathrm{ps}`$. A Kohlrausch-Williams-Watt fit ($`f(t)=\mathrm{exp}((t/A)^B)`$) from $`0`$ to $`80\mathrm{ps}`$ reveals a total correlation time of $`80\mathrm{ps}`$ by analytical integration . Contrary, the CMC II gyration radius decay has four domains: within the first picosecond there is a very fast and sharp decay (not resolved in figure 5) followed by a decay comparable to the first one for CMC I with a time constant of $`37\mathrm{ps}`$. The mid-time relaxation is in the same order of magnitude as for CMC I, but there is a vast difference in the long time decay, where CMC II has a very slow mode with a time constant of $`2200\mathrm{ps}`$. Averaging the whole range from $`0`$ to $`80\mathrm{ps}`$ with a Kohlrausch-Williams-Watt fit we calculate a total correlation time of $`900\mathrm{ps}`$.

The decay of the PAA gyration radius is less complex. It is almost monoexponential with a correlation time of $`110\mathrm{ps}`$. Compared with CMC – where monomers are larger and intramolecular hydrogen bonding is an issue (see section 3.3) – PAA lacks any fast decay. One may speculate about the reason for the different relaxation behaviour. Poly(acrylic acid) behaves monoexponentially as expected for a generic polymer chain in solution. The relaxation time of $`100\mathrm{ps}`$ is also found for both CMCs at some time scales which supports the idea that it is somehow characteristic for a polyelectrolyte of this size in water. The extremely slow long time relaxation (at least 2$`\mathrm{ns}`$) of CMC II is also easily explained by its rigid globular conformation. The origin of the short-time ($`t15\mathrm{ps}`$) relaxation behaviour of both CMCs (characteristic time $`40\mathrm{ps}`$), which is faster than for PAA, is as yet not clear. It is possibly due to the weaker hydration of CMC (chapter 3.3.2) which could decrease the short-time local friction experienced by the polymer.

### 3.2 Local Chain Properties

Two selected CMC I glucosidic torsional angle’s $`\varphi `$ (C4-O4-C1’-C2’) and $`\psi `$ (C3-C4-O4-C1’) trajectories are shown in figure 6. The positions of the maxima differ (for all angles) from results found for unsubstituted cellobiose in vacuum  and in water : For aqueous solutions there are maxima for ($`\varphi `$, $`\psi `$) at ($`284`$, $`63`$), ($`274`$, $`193`$), ($`294`$, $`242`$) and ($`50`$, $`237`$) degrees (in our definition). We found most of the linking dihedral angles to be unimodally distributed. Deviation from unimodality depend on the position along the chain and the substitution pattern: Torsion $`\psi `$ is unimodal for chain-terminating glucosidic links. A second $`\psi `$-state is populated only for links connecting repeat unit three with its neighbors. Occupancies of the individual states for CMC I are given in table 5. Since there are no significant barriers between the states, transitions or, rather, oscillations between them occur on a picosecond scale for all torsions. The terminal torsion $`\varphi `$ relaxes faster than the others ($`2\mathrm{ps}`$ rather than $`3`$$`4\mathrm{ps}`$). Other than this, no significant trends can be found. The relaxation time is calculated from the torsional angle auto correlation function $`\mathrm{\Phi }(t)\mathrm{\Phi }(0)`$ and a subsequent fit and integration.

There is hardly any torsional dynamics for CMC II, which is a consequence of the static folded state, with torsional links less mobile than those of CMC I. Example torsions are shown in figure 7.

Torsions are more complex for PAA (figure 8), with up to three states for each torsion. The time in which the torsions change ($`2`$$`5\mathrm{ps}`$) is in the same order of magnitude as for CMC I. There is no apparent correlation between the tacticity of adjacent PAA monomers and their backbone torsional statistics.

### 3.3 Hydrogen Bonding

For biopolymers and modified biopolymers hydrogen bonding (H-bonds) is one of the most significant interaction on atomistic scales, both intramolecular (conformation) and intermolecular (solvation, aggregation). Both CMC and PAA undergo intensive hydrogen bonding. In CMC hydroxy groups, carboxylic oxygens, and to a smaller extent, ether oxygens are involved. Poly(acrylic acid) can only form H-acceptor bonds via carboxylic groups.

Our definition of a hydrogen bond is geometric: If the two oxygens are within a distance of less than $`0.29\mathrm{nm}`$ with the hydrogen between them, so that the angle O-H-O is greater $`130^{}`$, a hydrogen bond is assumed. In our polyelectrolyte-water systems, we find both solute-solute, solute-solvent and solvent-solvent hydrogen bonding. The latter is only analyzed in the solute’s first hydration shell as defined from radial distribution functions (below).

The lifetimes $`\tau `$ of hydrogen bonds may be calculated from bond-correlation functions : For each possible hydrogen-bond we calculate an array, which contains ’true=$`1`$’ in case of a bond, otherwise ’false=$`0`$’. The time-correlation function is calculated on this array, integration leads to the lifetime $`\tau `$. The integration was done by fitting a stretched exponential (Kohlrausch-Williams-Watt) $`\mathrm{exp}[(t/\alpha )^\beta ]`$ in the range $`0`$ to $`80\mathrm{ps}`$ and followed analytical integration.

#### 3.3.1 Solute-Solute Hydrogen Bonding

The average absolute numbers of intramolecular hydrogen bonds are $`1.8`$ for CMC I and $`7.5`$ for CMC II respectively (see table 6 and figure 9). The large number of intramolecular hydrogen bonds stabilizes the much more compact structure of CMC II compared to CMC I.

For CMC, the hydrogen bonds lifetime $`\tau `$ differs considerably between the CMC I and CMC II molecule: CMC I has a bond lifetime of $`35\mathrm{ps}`$, CMC II of $`205\mathrm{ps}`$ (table 6), indicating that the cyclic compact structure of CMC II is stabilized by long-lived hydrogen bonds.

Intramolecular hydrogen bonds can be further separated into two subsets: *(a)* intra-ring and *(b)* inter-ring bonds, which connect two distinct glucose moieties and which are presumed to influence the overall chain structure of CMC: The ratio of inter- and intra-ring bonds for CMC I is $`3:1`$. Lifetimes are shorter for intra-ring bonds ($`\tau =23\mathrm{ps}`$), which demonstrates a lower stability, compared to inter-ring bond with $`\tau =40\mathrm{ps}`$. Further differences can be spotted in the hydrogen bonding atom types: intra-ring bonds involve O3-O6 hydrogen bonds ($`20\%`$) and to a larger amount O3-COO- bonds ($`70\%`$). Bonds between rings involve mainly O2-COO- ($`60\%`$) and O3-COO- and O6-COO- ($`15\%`$ each) oxygen atoms.

Thus the picture of intramolecular hydrogen bonds in CMC I is as follows: There is a high fraction of long-lived, ring-connecting hydrogen bonds, intra-ring bonds are of less importance. Inter-ring hydrogen bonding happens mostly via COO- groups, which makes *(a)* the interaction energy favourable and *(b)* involves a flexible, exocyclic group into binding. For this reason, changes in the global chain conformation may be compensated by carboxy methyl torsional changes to keep the OH-O alignment in an energetically low geometry. Intra-ring H-bonds, on the other hand, are either formed by neutral OH groups, which have a lower energy, or they necessitate sharp turns of the carboxy methyl group, which leads to distortions.

In CMC I, inter-ring hydrogen bonds are limited to neighboring rings (table 7). This is different for CMC II, where also more distant rings are connected by inter-ring hydrogen bonds: In contrast to CMC I, $`90\%`$ of the inter-ring hydrogen bonds connect rings which are not nearest neighbors: on average bonds span $`3`$$`4`$ rings. These bonds are found to be very stable, some of them are intermittently present for an overall time of $`1.7\mathrm{ns}`$, thus living for more than half of the simulation.

Furthermore, in CMC II there are several pairs of hydrogen bonds (about $`15`$), which bind almost exclusively (see figure 10). Out of these pairs, there is solely one bond present, but alternating with the second one. The presence of one bonds prevents the second bond. Even, if the mean lifetime of an individual bond is short, both bonds together form a steady connection between two glucose rings. One example is the bond pair O6(ring 3)-COO(ring 7) and O6(ring 5)-COO(ring 7), where the carboxylate oxygen flips between the two alcohol oxygens. Each bond lives $`1\mathrm{ns}`$ and together they span not less than $`2\mathrm{ns}`$. However, for $`1.7\mathrm{ns}`$ only one of them is exclusively present. Thus, for most of the time, one bond excludes the other. This evident for H-bonds formed via a COO- group, which account for every second ($`56\%`$) of all solute-solute bonds.

A distribution of all bond-bond cross correlation coefficients for CMC II is shown in figure 12. They are calculated as

$$c_{ij}=(h_i(t)h_i)(h_j(t)h_j)(h_i(t)h_i^2h_j(t)h_j^2)^{1/2}$$

where $`h_i(t)`$ is $`1`$ when hydrogen bond $`i`$ is present and $`0`$ otherwise. Most hydrogen bonds do not have a significant correlation: $`10`$ hydrogen-bond-pairs have a strong negative ($`c_{ij}<0.30`$), $`45`$ bond-pairs a positive correlation ($`c_{ij}>0.3`$). The anti-correlated bond pairs often involve one common oxygen of a COO--group, which alternates between two different hydrogen-donor groups.

High correlations ($`c_{ij}>0.6`$) occur particularly, if the compact configuration of the CMC II oligomer allows for some bonds to form simultaneously, like H-bonds involving O2H and O6H donor groups of one, and a single carboxylic group of a second ring. The COO- group is embedded between the two H-donors (figure 11). The effect is very strong for the flexible COO- group at rings 7 and 8, but is still noticable for other acceptors. Examples of hydrogen bond trajectories (Fig 13) illustrate the cases of (anti) correlated H-bonds.

The cross correlation coefficient $`c`$

$$c=\frac{(nn)(rR_{\mathrm{ete}})}{\sqrt{(nn)^2(R_{\mathrm{ete}}R_{\mathrm{ete}})^2}}.$$
(1)
shows some coincidence between the number of solute-solute hydrogen-bonds $`n`$ and the polyelectrolyte’s end-to-end distance $`R_{\mathrm{ete}}`$ for CMC I. Whereas there is hardly any correlation of extension and the total number of intra-CMC hydrogen-bonds ($`c=0.05`$), it is slightly more pronounced ($`c=0.16`$) for the number of intra-ring hydrogen bonds versus the end-to-end distance: If the molecular extension shrinks, the number of intramolecular hydrogen bonds increases. The compact polymer structure favours more contacts. Similar correlations for CMC II would not be very meaningful, because of the rigid structure which fixes the end-to-end distance.

#### 3.3.2 Hydrogen Bonds to Solvent

A general picture of the solvent distribution around CMC is given by radial distribution functions (RDF) between the center of mass (defined as in section 2) of the glucose unit and water oxygens (figure 14).

The differences in the peak heights can be described and explained as follows: The outermost glucose rings exhibit the best defined hydration shell. Environments of inner rings are more perturbed, either by neighboring rings or by attached carboxylic groups. This is evident for the doubly substituted ring 3 in CMC I, whose first solvation peak disappears entirely. The CMC II radial distributions for rings 1 and 3 show multiple peaks with the first maximum at approximately $`0.45\mathrm{nm}`$. Only ring 7 shows two different peaks centered at $`0.3\mathrm{nm}`$ and $`0.6\mathrm{nm}`$, where the first peak at $`0.3\mathrm{nm}`$ is caused by the solvation of the C6 carboxy methyl group. Thus from the appearance of the RDFs, the CMC I solvation shell seems to be more intact than that at CMC II. This is the consequence of the open-chain configuration of CMC I as opposed to the compact structure of CMC II, which excludes solvent molecules from the vicinity of different monomers to a different extent.

As PAA has a simple configuration, the solvation shell is more uniform: the chain ends attract more water than inner segments, but, apart from that, the center of mass-water RDF is the same with the first peak at $`0.35\mathrm{nm}`$.

On average the CMC I oligomer has $`13`$ donor and $`38`$ acceptor hydrogen bonds with water, CMC II $`9`$ donor and $`35`$ acceptor bonds (table 8). The difference is in line with the result for intramolecular bonds: The CMC II molecule, which heavily bonds with itself, has fewer free sites for hydrophilic interaction with the solvent. Poly(acrylic acid), with no hydrogen bond donors of its own, has on average no less than $`136`$ bonds with water. As an estimate, this corresponds to water molecules with a total mass of $`2300`$ amu, which is $`1.5`$ times the polymer’s own mass. The mass of associated water is smaller for CMC: About $`60\%`$ of the polymer’s mass for CMC I and only $`50\%`$ for CMC II. The detailed distributions of H-Bonds are shown in figure 15. They are approximately Gaussian in all cases and they reflect the trends already seen in the averages.

In order to compare the hydration of the different molecules on an equal basis, one may normalize the number of H-bonds by the mass of the polymer. This yields the H-bonds per weight, which should in some way be correlated with the specific enthalpy of solvation. One finds $`0.036`$ bonds/amu for CMC I, $`0.029`$ bonds/amu for CMC II and more than twice as many for PAA: $`0.083`$ bonds/amu. Thus there seems to be a competition between inter- and intra-molecular hydrogen bonds for carboxy methyl cellulose, which is dominated by the intramolecular bonds in case of CMC II and to a lesser extent for CMC I. In PAA, there is no competition from intramolecular hydrogen-bonds. In addition, one can reasonably expect that hydrogen-bonds to COO- with its negative charge are energetically more favorable than to –OH or –O–. Hence, PAA should have a more negative enthalpy of solvation than CMC. We are, however, only aware of experimental data for PAA to compare with . Klein et al. used solution calorimetric measurements to determine the solution enthalphy of about $`25\mathrm{kJ}`$/(mol repeat unit) ($`0.27\mathrm{kJ}/\mathrm{g}`$) for sodium-PAA, this value corresponds to $`575\mathrm{kJ}/\mathrm{mol}`$ for the PAA 23-mer.

### 3.4 Structure of the solvation shell

The arrangement of water molecules in the vicinity of the polymer chains has been investigated. To this end, all water molecules within $`0.6\mathrm{nm}`$ of any polymer atom were sorted into clusters. A water molecule was taken to be part of a cluster if it formed a hydrogen bond to any other water molecule of that cluster. Hence, the assignment of molecules to clusters depends on our geometrical definition of a hydrogen bond (Sect. 3.3).

For the three polymers, examples of water clusters in the immediate vicinity of the solute are shown in Figure 16. For both CMCs, such clusters consist of $`2`$$`5`$ molecules and are attached to the solute via $`2`$$`3`$ hydrogen bonds (in addition to individual water molecules hydrogen-bonded to CMC which are not shown). Depending on their size, the clusters can connect two (CMC I: a, b; CMC II: a-c) or more (CMC I: c,d) cellulose rings. Poly(acrylate) is quite different due to the shorter distances between repeat units. Typically, one finds single water molecules, that are not part of water clusters and that bridge neighboring carboxylate groups (Fig 16, PAA I). The clusters are, however, fragile arrangements undergoing rapid exchange (order of picoseconds, Fig 17) with other water molecules and they do provide little if any stabilisation to particular CMC conformations. Only a small difference is found between both CMC oligomers (figure 17), whereas clusters in the hydration sphere of PAA tend to live longer. This could be due to the high concentration of negative charges on the PAA backbone, which slow down the general motion of hydration water by their strong electrostatic field.

### 3.5 Counterions

The distribution of counterions around the solute molecules has also been investigated. The solutes are fragments of polyelectrolytes, and since the work of Manning  the distribution of counterions in the field of polyelectrolytes has been of considerable theoretical interest (see, e.g., ref.  and references therein). Because of the small size of the charged solutes investigated here and the small number of counterions, this study has to be limited to the molecular neighborhood of the solutes.

Radial distribution functions (RDFs) between the centres of mass of the repeat units (for definition, see Sect. 2) and the Na+ are summarised in Figure 18 a. Both CMCs are similar in shape, but different in intensity: They have a close peak at about $`0.4\mathrm{nm}`$ and a broad second peak centered at about $`0.7\mathrm{nm}`$. Figure 18 b shows that the first peak ($`r<0.5\mathrm{nm}`$) is dominated by Na+ approaching the O2 oxygen of the cellulose moiety, whereas the second peak ($`0.5`$$`1\mathrm{nm}`$) arises from Na+ ions near the carboxylic groups which, in CMC, are more remote from the centre of mass. Both the (dynamic) flexibility of the carboxylate and the (static) different substitution patterns on different AGUs account for the large width of the second peak. CMC II has a lower intensity, especially at short distance. This arises from its more compact structure which prevents counterions from approaching it freely from all sides. This behaviour is similar to what was already seen for the solute-solvent H-bonds.

The RDF of poly(acrylate) is better defined due to the larger number of counterion-monomer pairs. It shows a relatively sharp first maximum at about $`0.5\mathrm{nm}`$ caused by Na+ ions near the carboxylate and a broader second maximum ($`0.8\mathrm{nm}`$) possibly due to solvent-separated ion pairs.

It is interesting to note that on average the CMC oligomer keeps much fewer counterions ($`0.41`$) in its electrostatic vicinity (distance to any atom of the polymer lower than the electrostatic screening length of the solvent $`<0.64\mathrm{nm}`$) than the PAA oligomer ($`13`$), in spite of both being of similar molecular weight. This is caused by the difference in charge.

## 4 Conclusions

Even though poly(acrylate) and carboxy methyl cellulose both are water-soluble polyelectrolytes, their behaviour in water and towards water differs markedly. This is due to the different charge density as well as to the different type and quality of hydrogen bonds that either forms with water. In PAA, there is one strong hydrogen bond with the deprotonated carboxylate acting as an acceptor. In CMC, the smaller density of carboxylates is only partly set off by the possibility of forming both donor and acceptor hydrogen bonds to the alcoholic OH groups of the cellulose. (Hydrogen bonds to the ether oxygens are irrelevant.) Taken per molecular weight of the polymer, it seems safe to say that PAA forms at least twice as many hydrogen bonds as CMC and that they are of larger binding energy (charge-dipole, rather than, dipole-dipole). Based on this argument, the solvation of PAA in water should be more exothermic than that of CMC. Unfortunately, no measurements appear to be available for comparison.

The comparison of the two CMC oligomers shows that the particular carboxy methylation pattern has an immense influence on the local structure in solution. The two assume entirely different conformations: CMC I is stretched and flexible, whereas CMC II remains in a rigid cyclic conformation. While it cannot be ruled out completely that one molecule may not have reached equilibrium and may still be stuck in a local minimum, this is unlikely in view of the long simulation times of several nanoseconds. If there are metastable conformations they have to be long-lived. We are therefore left to conclude that industrial CMC with its statistical substitution of OH groups, behaves locally very diversely. As a consequence of its globular structure, CMC II shows more intramolecular hydrogen bonds than CMC I, fewer hydrogen bonds to water, slower hydrogen bond dynamics, and fewer contacts with the counterions.

## 5 Acknowledgements

We would like to thank the John von Neumann Computer Center of the Forschungszentrum Jülich for Cray T90 time. Fruitful discussions with Markus Deserno and Roland Faller are gratefully acknowledged.