text:
# Self-Consistent Theory of Rupture by Progressive Diffuse Damage

## 1 Introduction

The fracture of materials is a catastrophic phenomenon of considerable technological and scientific importance. Despite the large amount of experimental data and the considerable effort that has been undertaken by material scientists , many questions about fracture remain standing. There is no comprehensive understanding of rupture phenomena but only a partial classification in restricted and relatively simple situations. This lack of fundamental understanding is reflected in the absence of reliable prediction methods for rupture based on a suitable monitoring of the stressed system.

Some progresses have been obtained in recent years in the Physics community. Based on analogies with phase transitions, several groups - have proposed that, in heterogeneous materials with disorder such as fiber composites, rocks, concrete under compression and materials with large distributed residual stresses, rupture is a genuine critical point, i.e., the culmination of a self-organization of diffuse damage and micro-cracking characterized by power law signatures. Experiments , numerical simulations  and theory  confirm this concept.

As a signature of criticality, acoutic emissions radiated during loading exhibit an acceleration of their rate close to rupture . Specifically, under a constant stress rate, the cumulative acoustic energy $`E(t)`$ released up to time $`t`$ can be expressed as

$$E(t)=E_0B(t_ct)^\alpha ,$$
(1)
with $`B>0`$ and $`0<\alpha <1`$. Expression (1) corresponds to a rate $`dE/dt`$ of acoustic energy release diverging at the critical rupture time $`t_c`$. This behavior (1) has been at the basis of previous claims that rupture is a critical phenomenon. In addition, this power law (1) as well as extensions with log-periodic corrections have been suggested to be useful for prediction .

Our purpose here is to present, extend and analyze a simple self-consistent model of damage that predicts a behavior similar to (1). We explore its different regimes and then improve on its “mean-field” version which predicts an unrealistic finite-time singularity. In this goal, we propose to use the general functional renormalization approach developed by Yukalov and Gluzman to cure this anomaly. We show how this technique allows us to change an unrealistic singularity into the observed behavior (1) with a reasonable exponent $`z=1/2`$, without introduction of any extra parameters in the theory.

## 2 Cumulative damage model

Initially introduced as a global “mean field” (uniform) description of the global deterioration of the system at the macroscopic scale , the concept of “damage” has been extended at the mesoscopic scale to describe the heterogeneity and spatial variability of damage in different locations within the material . We use the formulation of Zobnin  and Rabotnov  to show how it leads naturally to a finite-time singularity. We first recall briefly the integral formulation of Rabotnov  (pages 166-170) and then transform it in differential form to exhibit the fundamentally nonlinear geometrical origin of the singularity.

A material is subjected to a stress $`\sigma _0`$ at large scale and each point $`r`$ within it carries a damage variable $`d(r,t)`$. When $`d`$ reaches the threshold $`d^{}`$ as some location, this local domain is no more able to sustain stress and a microcrack appears, leading to a redistribution of the stress field around it according to the laws of elasticity. The local damage $`d(r,t)`$ at point $`r`$ at time $`t`$ is supposed to evolve in time according to

$$\frac{d(d)}{dt}=[\sigma (r,t)]^m,$$
(2)
where $`\sigma (r,t)`$ is the local stress field at point $`r`$ at time $`t`$ and $`m`$ is a damage exponent which can span values from $`0`$ to close to $`\mathrm{}`$ depending upon the material. In the discrete 2D models of Refs., it was shown that rupture reduces to the percolation model in the limit $`m0`$. In the other limit, $`m\mathrm{}`$, rupture occurs through a one crack mechanism.

Following Rabotnov , we assume that a major crack dominates the rupture process. If only one crack is present within the system, the stress $`\sigma (r,t)`$ is easily calculated. Considering only the possibility of a linear straight crack of half-length $`a(t)`$ advancing within the material at velocity $`da/dt`$ (see  for generalizations to self-affine crack geometries), it is enough to calculate the stress field on the points ahead of the crack to fully characterize the rupture dynamics. For a planar elastic material subjected to a uniformly distributed antiplane stress at infinity with a crack lying on the $`y`$-axis between $`a(t)`$ and $`a(t)`$, the stress field at point $`y`$ on the $`y`$-axis beyond the crack tip is

$$\sigma (z,t)=\frac{2\sigma _0}{3}\frac{z}{\sqrt{z^2[a(t)]^2}}.$$
(3)
The mean-field approximation made in this first version of the model consists in assuming that the stress field is not modified by the non-vanishing and non-uniform damage field. This means that the elastic coefficients are taken constant and independent of the progressive damage, except of course when the damage reaches its rupture threshold $`d^{}`$.

The law describing the growth of the crack, i.e., the dynamics $`a(t)`$, is obtained from the self-consistent condition that the time it takes from a point at $`y`$, at the distance $`ya(\tau )`$ from the crack tip at time $`\tau `$, for its damage to reach the rupture threshold $`d^{}`$ is exactly equal to the time taken for the crack to grow from size $`a(\tau )`$ to the size $`a(t)=y`$ so that its tip reaches the point $`y`$ exactly when it ruptures. This is illustrated in figure 1. Mathematically, this self-consistent condition is that the integral of (2) from time $`0`$ at which the pre-existing damage was $`0`$ till time $`t`$ at which the crack tip passes through $`y`$ is such that $`d`$ reaches exactly the threshold $`d^{}`$ at the time $`t`$. Two conditions must thus be verified simultaneously:

1. $`a(t)=y`$ (the crack tip reaches point $`y`$) and
2. $`d(y,t)=d^{}`$ (the damage at $`y`$ reaches the rupture threshold).

## 3 The linear damage law: $`m=1`$

We first consider the linear damage law $`m=1`$ corresponding to the initial formulation of Zobnin . This case has also been investigated and solved in  in the context of crack growth due to electromigration rather than mechanical stress (the current plays the role of the stress and, in the antiplane case studied here, the two problems are formally identical). This model is particularly interesting since it allows both for an exact solution and an exact renormalization in the functional renormalization scheme . It also provides a benchmark for approximate solutions in the general case $`0<m<2`$ as we discuss below.

We now proceed to give the equation for the crack dynamics and its solution. By integration of (2), the two self-consistent conditions expressed for the case $`m=1`$ lead to

$$_0^t𝑑\tau \frac{2\sigma _0}{3}\frac{a(t)}{\sqrt{[a(t)]^2[a(\tau )]^2}}=d^{},$$
(4)
where the loading stress $`\sigma _0`$ can depend on time. The solution of this integral equation provides the time evolution $`a(t)`$ of the macro-crack. To get it explicitely, we set

$$z=[a(t)]^2\mathrm{and}\zeta =[a(\tau )]^2.$$
(5)
Changing the variable of integration from $`\tau `$ to $`\zeta `$ gives

$$_{z_0}^z𝑑\zeta \sigma _0\frac{(d\tau /d\zeta )}{\sqrt{z\zeta }}=\frac{3}{2}\frac{d^{}}{\sqrt{z}}.$$
(6)
This equation (6) is an Abel equation with index $`1/2`$, involving a fractional integral operator . Defining the Abel operator $`I_\alpha ^{}`$ acting on the function $`f(t)`$ as

$$I_\alpha ^{}\{f\}=_0^t\frac{(ts)^\alpha }{\mathrm{\Gamma }(1+\alpha )}f(s)𝑑s,$$
(7)
the product of two such Abel operators is

$$I_\alpha ^{}I_\beta ^{}\{f\}=I_\alpha ^{}\{I_\beta ^{}\{f\}\}=_0^t𝑑\tau \frac{(t\tau )^\alpha }{\mathrm{\Gamma }(1+\alpha )}_0^\tau \frac{(\tau s)^\beta }{\mathrm{\Gamma }(1+\beta )}f(s)𝑑s$$

$$=_0^t𝑑sf(s)\frac{1}{\mathrm{\Gamma }(1+\alpha )\mathrm{\Gamma }(1+\beta )}_s^t𝑑\tau (t\tau )^\alpha (\tau s)^\beta =_0^t𝑑sf(s)\frac{(ts)^{\alpha +\beta +1}}{\mathrm{\Gamma }(2+\alpha +\beta )}.$$
(8)
This shows that

$$I_\alpha ^{}I_\beta ^{}\{f\}=I_{\alpha +\beta +1}^{}\{f\}.$$
(9)
We thus see that $`I_\alpha ^{}I_{1\alpha }^{}\{f\}=I_0^{}`$, which is nothing but the integral operator. The inverse of the Abel operator $`I_\alpha ^{}`$ is thus $`\frac{d}{dt}\left(I_{1\alpha }^{}\right)`$. Applying this result to (6), we find

$$\frac{d\tau }{dz}=\frac{3}{2}\frac{d^{}}{\pi \sigma _0}\frac{d}{dz}_{z_0}^z\frac{d\zeta }{\sqrt{\zeta (z\zeta )}}.$$
(10)
Calculating the integral in the r.h.s. of (10), performing the derivative and inverting to get $`dz/d\tau `$, we get

$$\frac{dz}{d\tau }=\frac{2\pi \sigma _0}{3d^{}}\frac{z\sqrt{zz_0}}{\sqrt{z_0}}.$$
(11)
Replacing $`z`$ by $`[a(t)]^2`$ leads to the differential equation for the crack half-length $`a(t)`$

$$\frac{da}{dt}=\frac{\pi \sigma _0}{3d^{}}a\sqrt{\left(\frac{a}{a_0}\right)^21},$$
(12)
which is exactly equivalent to the self-consistent integral equation (4). It is remarkable that the local growth equation (12) embodies exactly the same physics as the long-term memory integral (4).

For simplicity, let us take the loading stress $`\sigma _0`$ constant. This situation is generic of experiments measuring the lifetime of structures under a constant load. At sufficiently long times for which $`a(t)>>a_0`$, expression (12) reduces to

$$\frac{da}{dt}\frac{\pi \sigma _0}{3a_0d^{}}a^2.$$
(13)
Equation (13) is characteristic of a solution going to infinity in finite time. Indeed, we can write (13) as $`da/dtra`$, with a growth rate $`ra`$. The generic consequence of a power law acceleration in the growth rate $`ra^\delta `$ with $`\delta >0`$ is the appearance of a singularity in finite time:

$$a(t)(t_ct)^\beta ,\mathrm{with}\beta =\frac{1}{\delta }\text{ and }t\text{ close to }t_c.$$
(14)
Equation (13) is said to have a “spontaneous” or “movable” singularity at the critical time $`t_c`$ , the critical time $`t_c`$ being determined by the constant of integration, i.e., the initial condition $`a(t=0)=a_0`$. Note the intriguing fact that the $`(t_ct)^1`$ singularity appears as the solution of a linear mechanical problem. The source of the quadratic nonlinearity is the non-local geometrical condition that the delayed action of the stress field on the cumulative damage should coincide exactly with the passage of the crack tip. The nonlinear finite-time singularity has thus fundamentally a non-local geometrical origin, or alternatively can be seen to result from a long-term memory effect.

The exact solution of (12) is easily obtained by integration:

$$a(t)=\frac{a_0}{\mathrm{cos}\left(\frac{\pi \sigma _0}{3d^{}}t\right)}.$$
(15)
This retrieves the solution obtained by Zobnin  and Rabotnov . We verify directly that the singularity occurs when the cosine goes to zero, i.e., when the argument reaches $`\pi /2`$, i.e., for $`t_c=3d^{}/2\sigma _0`$. Since the cosine vanishes linearly with time, this recovers the asymptotics (14) with the exponent $`\beta =1`$, as predicted by the asymptotic equation (13).

## 4 The nonlinear damage law with finite-time singularity: $`0<m<2`$

### 4.1 Derivation of the differential equation for the crack dynamics

The case where $`0<m<2`$ can be similarly treated and our results here extend those of Zobnin  and Rabotnov . Our results retrieve those found in , obtained in the context of crack growth due to electromigration. For completeness and coherence in notation, we briefly present the method and the results which are more focused on the finite-time singularity.

For simplicity, we impose $`\sigma _0`$ constant. Integrating (2) and applying the self-consistent conditions 1-2 leads to

$$_0^t𝑑\tau \left(\frac{2\sigma _0}{3}\right)^m\frac{[a(t)]^m}{\left([a(t)]^2[a(\tau )]^2\right)^{\frac{m}{2}}}=d^{}.$$
(16)
We set again the change of variables (5) and changing the variable of integration from $`\tau `$ to $`\zeta `$ gives

$$_{z_0}^z𝑑\zeta \frac{(d\tau /d\zeta )}{\left[z\zeta \right]^{\frac{m}{2}}}=\left(\frac{3}{2\sigma _0}\right)^m\frac{d^{}}{z^{\frac{m}{2}}}.$$
(17)
This equation (6) is again an Abel equation with index $`m/2`$ if $`0<m<2`$.

In order to transform it into differential form, we could use the formalism of Abel operators. We choose a more transparent and direct approach which is closely related. First, we multiply both sides of (17) by $`1/(yz)^{1\frac{m}{2}}`$ and integrate over $`z`$ from $`z_0`$ to $`y`$:

$$_{z_0}^y𝑑z\left(_{z_0}^z𝑑\zeta \frac{(d\tau /d\zeta )}{\left(z\zeta \right)^{\frac{m}{2}}(yz)^{1\frac{m}{2}}}\right)=\left(\frac{3}{2\sigma _0}\right)^md^{}_{z_0}^y\frac{dz}{z^{\frac{m}{2}}(yz)^{1\frac{m}{2}}}.$$
(18)
Changing the order of integration in the l.h.s. leads to

$$_{z_0}^z𝑑\zeta \frac{d\tau }{d\zeta }\left[_\zeta ^y\frac{dz}{\left(z\zeta \right)^{\frac{m}{2}}(yz)^{1\frac{m}{2}}}\right]=\left(\frac{3}{2\sigma _0}\right)^md^{}_{z_0}^y\frac{dz}{z^{\frac{m}{2}}(yz)^{1\frac{m}{2}}},$$
(19)
where we have used the equality of the triangle $`_{z_0}^y𝑑z_{z_0}^z𝑑\zeta =_{z_0}^z𝑑\zeta _\zeta ^y𝑑z`$.

The integral in the square bracket in the l.h.s. of (19) can be expressed through the Euler beta-function $`B(a,b)`$:

$$_\zeta ^y\frac{dz}{\left(z\zeta \right)^{\frac{m}{2}}(yz)^{1\frac{m}{2}}}=B(1\frac{m}{2},\frac{m}{2})=\frac{\mathrm{\Gamma }\left(1\frac{m}{2}\right)\mathrm{\Gamma }\left(\frac{m}{2}\right)}{\mathrm{\Gamma }\left(1\frac{m}{2}+\frac{m}{2}\right)}=\frac{\pi }{\mathrm{sin}\left(m\frac{\pi }{2}\right)}.$$
(20)
We thus obtain

$$\frac{\pi }{\mathrm{sin}\left(m\frac{\pi }{2}\right)}_{z_0}^z𝑑\zeta \frac{d\tau }{d\zeta }=\left(\frac{3}{2\sigma _0}\right)^md^{}_{z_0}^y\frac{dz}{z^{\frac{m}{2}}(yz)^{1\frac{m}{2}}}.$$
(21)
After differentiation with respect to $`z`$, we get

$$\frac{d\tau }{dz}=\left(\frac{3}{2\sigma _0}\right)^md^{}\frac{\mathrm{sin}\left(m\frac{\pi }{2}\right)}{\pi }\frac{d}{dz}_{z_0}^z\frac{d\zeta }{\zeta ^{\frac{m}{2}}(z\zeta )^{1\frac{m}{2}}}.$$
(22)
### 4.2 Asymptotic solution close to the finite-time singularity

The solution of (22) can be obtained for large crack sizes $`a(t)`$, i.e., large $`z`$. In this goal, we replace the term $`(z\zeta )^{1\frac{m}{2}}`$ in the integral in the r.h.s. of (22) by $`z^{1\frac{m}{2}}`$, neglecting $`\zeta `$ compared to $`z`$. Intuitively, this is justified over the whole domain of integration because the contribution from the domain where $`\zeta `$ is not negligible compared to $`z`$ is finite, since the power $`1\frac{m}{2}`$ is less than one, corresponding to an integrable singularity.

With this approximation, the integral can be performed, the derivative taken and after inverting, we get

$$\left(\frac{3}{2\sigma _0}\right)^md^{}\frac{\mathrm{sin}\left(m\frac{\pi }{2}\right)}{\pi }\frac{dz}{d\tau }=z_0^{1+\frac{m}{2}}z^{2\frac{m}{2}}.$$
(23)
Using $`a=\sqrt{z}`$ as defined in (5), we obtain

$$a(t)=\frac{a_0}{\left(1\frac{t}{t_c}\right)^\beta },$$
(24)
where

$$t_c=\left(\frac{3}{2\sigma _0}\right)^md^{}\frac{2\mathrm{sin}\left(m\frac{\pi }{2}\right)}{\pi (2m)},$$
(25)
and

$$\beta =\frac{1}{2m}.$$
(26)
Note that the exact asymptotics (14) of the case $`m=1`$ previously solved exactly is recovered, with the correct exponent $`z(m=1)=1`$ and a rather good approximation of the critical $`t_c`$: while the exact value is $`t_c=3d^{}/3\sigma _0`$, expression (25) predicts $`\frac{2}{\pi }t_c`$, i.e., $`36\%`$ lower. The critical time $`t_c`$ as a function of $`m`$ is smooth with no accident or divergence over the whole interval. In particular, the estimated critical time for the limit $`m2^{}`$ is equal to $`9d^{}/4\sigma _0^2`$.

In contrast, the exponent $`z`$ increases from $`z(m0^+)=1/2`$ to $`+\mathrm{}`$ as the damage exponent $`m`$ varies from $`0`$ to $`2`$. The limit $`z(m0^+)=1/2`$ can be rationalized as follows. This limit $`m0^+`$ corresponds to the situation where damage becomes independent of stress. As a consequence, reintroducing some heterogeneity for instance on the pre-existing damage, rupture is then equivalent to percolation, as the parts of the system that break as a function of time are determined by the damage accumulating at the same rate for all point but with different random initial values. In mean field percolation  obtained through the consideration of one-dimensional percolating paths consistent with the present one-crack geometry, the elastic energy under constant load diverges a $`(t_ct)^1`$ where $`1`$ is the mean field value of the exponent $`t`$ for conductivity (which is the same as elasticity in the scalar mode III version of mechanical deformations used here). Since the elastic energy is proportional to the square of the crack length, we get the prediction $`a(t)(t_ct)^{1/2}`$. This reasoning holds if the exponent is a smooth function of disorder and geometry (the present studied here is a the zero-disorder limit).

The divergence of $`z`$ at $`m=2`$ signals a change of regime that we study in the next section.

## 5 The nonlinear damage law with $`m2`$

For $`m2`$, the integrals in the equations (16) and (17) diverge at $`a(\tau )=a(t)`$, since the negative power with exponent $`m/2`$ is no more integrable. Technically, the main difference between the cases $`m<2`$ and $`m2`$ is that we need to regularize the infinity in the expression (20) by introducing some sort of dimensionless cut-off. The important physical message is that the regime where $`m2`$ is controlled by a novel physical parameter, which we identify as a length scale associated with the damage law. In other words, the physics of the rupture is inherently controlled by the choice of the cut-off, i.e., by the existence of a microscopic length scale. We could summarize the situation by saying that there is no continuous limit to the theory for $`m2`$. This is similar to previous observations obtained in a dynamical theory of rupture front propagation . We now present two ways for regularizing the divergence and thus for obtaining a meaningful theory of rupture.

### 5.1 Regularization by damage saturation at a microscopic scale

Before describing the physical content of the regularization we propose, we need to express the problem in a more manageable mathematical form. Since the culprit for the divergence is the integral (20) and the divergence occurs for $`\zeta z^{}`$, we introduce the variable

$$Z=\frac{z\zeta }{y\zeta },$$
(27)
and rewrite (20) as

$$_\zeta ^y\frac{dz}{\left(z\zeta \right)^{\frac{m}{2}}(yz)^{1\frac{m}{2}}}=_0^1𝑑ZZ^{\frac{m}{2}}(1Z)^{\frac{m}{2}1},$$
(28)
which makes apparent that the divergence is due to $`Z^{\frac{m}{2}}`$ at the lower bound $`0`$. It is thus natural to regularize by introducing a dimensionless cut-off $`ϵ>0`$ and replace (28) by

$$_ϵ^1𝑑ZZ^{\frac{m}{2}}(1Z)^{\frac{m}{2}1}b(m,ϵ).$$
(29)
The function $`b(m,ϵ)`$ is such that

$$\mathrm{lim}_{ϵ0^+}b(m,ϵ)=B(1\frac{m}{2},\frac{m}{2}),\mathrm{for}0<m<2,$$
(30)
where the beta function $`B(1\frac{m}{2},\frac{m}{2})`$ has been defined in (20).

In constrast, we have

$$b(m,ϵ)\frac{1}{ϵ^{\frac{m2}{2}}},\mathrm{for}m>2,$$
(31)
and

$$b(m,ϵ)\mathrm{ln}\frac{1}{ϵ},\mathrm{for}m=2,$$
(32)
showing that the divergence of the integral (28) is now encapsulated in the dependence of the factor $`b(m,ϵ)`$ on $`ϵ`$. This regularization scheme thus relies on the existence of the definite integral (29), by analogy to the case $`m<2`$.

Using the regularization (29), we obtain

$$\frac{d\tau }{dz}=\left(\frac{3}{2\sigma _0}\right)^m\frac{d^{}}{b(m,ϵ)}\frac{d}{dz}_{z_0}^z\frac{d\zeta }{\zeta ^{\frac{m}{2}}(y\zeta )^{1\frac{m}{2}}},$$
(33)
which extends (10) to the regime $`m2`$. Its formal solution obtained in implicit form is

$$t=\left(\frac{3}{2\sigma _0}\right)^m\frac{d^{}}{b(m,ϵ)}_{z_0}^z𝑑y\frac{d}{dy}_{z_0}^y\frac{d\zeta }{\zeta ^{\frac{m}{2}}(y\zeta )^{1\frac{m}{2}}}.$$
(34)
This regularization scheme allows to obtain exact solutions for integer $`m`$’s. We examine the solutions for $`m=2,3,4`$ and $`5`$ and then the general case. For m=2, we have the expression for all times given by

$$z_2(t)=z_0e^{t/t_{ϵ,2}},$$
(35)
where

$$t_{ϵ,2}=\frac{9d^{}}{4\sigma _0^2}\frac{1}{\mathrm{ln}\frac{1}{ϵ}}.$$
(36)
For $`m=3`$, $`z_3(t)`$ is the solution to

$$\sqrt{(z/z_0)1}\mathrm{tan}^1(\sqrt{(z/z_0)1})=\frac{t}{t_{ϵ,3}},$$
(37)
where

$$t_{ϵ,3}=\frac{27d^{}}{8\sigma ^3\sqrt{z_0}}ϵ^{\frac{1}{2}}.$$
(38)
For large times, we get

$$z_3(t)z_0\left(\frac{t}{t_{ϵ,3}}\right)^2.$$
(39)
For $`m=4`$, $`z_4(t)`$ is the solution to

$$\mathrm{ln}(z/z_0)+(z/z_0)1=\frac{t}{t_{ϵ,4}},$$
(40)
where

$$t_{ϵ,4}=\frac{81d^{}}{16\sigma ^4}ϵ.$$
(41)
For large times, we get

$$z_4(t)z_0\frac{t}{t_{ϵ,4}}.$$
(42)
For $`m=5`$, $`z_5(t)`$ is the solution to

$$\left((z/z_0)1\right)^{3/2}3\left((z/z_0)1\right)^{1/2}+3\mathrm{tan}^1\left[\left((z/z_0)1\right)^{1/2}\right]=\frac{t}{t_{ϵ,5}},$$
(43)
where

$$t_{ϵ,5}=\frac{243d^{}}{32\sigma ^5}ϵ^{3/2}.$$
(44)
For large times, we get

$$z_5(t)z_0\left(\frac{t}{t_{ϵ,5}}\right)^{2/3}.$$
(45)
More generally, at large times

$$z_m(t)z_0\left(\frac{t}{t_{ϵ,m}}\right)^{\frac{2}{m2}},$$
(46)
$$a_m(t)a_0\left(\frac{t}{t_{ϵ,m}}\right)^{\frac{1}{m2}},$$
(47)
where

$$t_{ϵ,m}ϵ^{\frac{m2}{2}}\frac{1}{b(m,ϵ)}.$$
(48)
From these solutions, it is apparent that the dynamics $`z(t)[a(t)]^2`$ is controlled by the characteristic time $`t_{ϵ,m}`$ defined in (48). Note that the inverse dependence of $`t_{ϵ,m}`$ on $`b(m,ϵ)`$ is obvious from the expression (34). As the cut-off $`ϵ0`$, $`t_{ϵ,m}0`$ and the global rupture occurs in vanishing time. The physical explanation of this phenomenon is as follows. For $`m2`$, the driving force $`\sigma ^m`$ of the damage law (2) is so strong close to and at the crack tip, that it takes effectively zero time for a point to be brought to the damage threshold. To see this, let us truncate the integral in (16) such that the upper bound is changed from $`t`$ to $`t\eta `$. The divergence of the integral at the crack tip means that the contribution to the cumulative damage occurring in the time interval from $`t\eta `$ to $`t`$ is larger (actually infinitely larger) than the contribution from time $`0`$ to time $`t\eta `$. This means that the progressive damage leading to the acceleration of $`a(t)`$ for $`m<2`$ is replaced by an infinite velocity as soon as we start from a finite crack and do not introduce the finite cut-off length.

This clarifies the physical meaning of the cut-off $`ϵ`$ defined in (29). A non-zero $`ϵ`$ means that the integral over $`\zeta `$ in (28) does not go all the way up to $`z`$. Translated in terms of physical distances, it means that the integral in $`\zeta `$ in the l.h.s. of (17) also does not go all the way up to $`z`$. Physically, this means that the damage on a given point ahead of the crack tip reaches the critical value $`d^{}`$ before the crack tip reaches that point. The value $`d^{}`$ is no more the rupture threshold but a saturation value. The crack tip dynamics is now determined by the condition that the damage at any given point $`y`$ reaches this saturation value $`d^{}`$ when the crack tip is at a fixed distance $`ϵ`$ from $`y`$. This condition embodies the existence of a microscopic length scale $`ϵ`$ such that the damage is no more defined as smaller scales.

The main result of our analysis is that the characteristic time scale $`t_{ϵ,m}`$ of the crack dynamics is controlled by the microscopic length scale. The theory has thus fundamentally no continuous limit. It is one of several interesting and important examples in physics where the macroscopic physics is completely controlled by the microscopic physics (the ultra-violet cut-off). This situation is found in many physical problems, for instance in correlation functions in two-dimensional systems , in non-linear diffusion  as well as in quantum electrodynamics . Note however the difference between the last two examples and the former ones: in our rupture problem as well as in the case of correlation functions in 2D systems, the ultraviolet cut-off appears naturally, as an atomic distance, while in the last two case, there is no meaningful natural cut-off, hence necessity to “cover-up” divergencies by the “renormalization” procedure .

### 5.2 Regularization by stress saturation at a microscopic scale

The previous regularization scheme invokes a saturation of the damage at a microscopic length $`ϵ`$. Alternatively, the saturation can occur on the stress field, whose mathematical divergence is bound to be rounded off at atomic scales. This provides another regularization scheme. To implement it, we use the continuous expression (3) for all distances from the crack tip down to a regularization length $`\mathrm{}`$ such that, for distances from the crack tip from $`0`$ to $`\mathrm{}`$, the stress is constant equal to $`\sigma (\mathrm{})`$ given by (3) with $`y=a(t)+\mathrm{}`$. This regularization is standard in the theory of damage and of plasticity. The idea is that a sufficiently large damage exerts a feedback on the stress field which then departs from its damage-free continuous expression (3). This extension to Rabotnov’s treatment provides a natural way for constructing a self-consistent theory of damage: not only does rupture occur by the cumulative effect of damage, damage has also the effect of smoothing out the mathematical singularity at the crack tip. The cut off $`\mathrm{}`$ has the physical meaning of a so-called process zone or damage zone and its introduction is fully consistent with the dynamical damage law (2). However as we have seen above, for a damage law with an exponent $`0<m<2`$, the resulting dynamics becomes insensitive to the existence of a microscopic length scale in the limit where it is small. In this sense, the regime $`0<m<2`$ is more universal and has a continuous limit.

We propose two models that implement these ideas.

#### 5.2.1 Saturation of the stress at a fixed distance to the crack tip

The regularization scheme used here is such that the stress is assumed to saturate at a fixed distance $`\mathrm{}`$ from the crack tip. Thus, for times $`\tau `$ up to $`t_{\mathrm{}}(t)`$, the damage at a fixed point that will be reached by the crack tip at time $`t`$ is growing under the influence of the stress field created by the crack. From time $`t_{\mathrm{}}(t)`$ up to time $`t`$, the damage is increasing linearly with time, since the stress is assumed constant and equal to the value it reaches at time $`t_{\mathrm{}}(t)`$. The saturation time $`t_{\mathrm{}}(t)`$ is determined by the equation

$$a(t)a(t_{\mathrm{}})=\mathrm{}.$$
(49)
In this version of the regularized theory, expression (16) is changed into

$$_0^t_{\mathrm{}}𝑑\tau \left(\frac{2\sigma _0}{3}\right)^m\frac{[a(t)]^m}{\left([a(t)]^2[a(\tau )]^2\right)^{\frac{m}{2}}}+\left(\frac{2\sigma _0}{3}\right)^m\frac{[a(t)]^m}{\left[(2a(t)\mathrm{})\mathrm{}\right]^{\frac{m}{2}}}(tt_{\mathrm{}})=d^{},$$
(50)
where $`t_{\mathrm{}}`$ is the time at which the crack tip is at the distance $`\mathrm{}`$ from the position it will have at time $`t`$ (see (49)).

Note that we now have two equations for two unknown $`a(t)`$ and $`a(t_{\mathrm{}})`$. The second term in the l.h.s. of (50) expresses the linear increase in damage from $`t_{\mathrm{}}`$ to $`t`$ under the saturated stress. In this model, $`\mathrm{}`$ is fixed and $`t_{\mathrm{}}`$ ajusts itself. The value of the saturated stress is not a constant but increases as the crack gets larger and larger, since it corresponds to the value at a fixed distance $`\mathrm{}`$ from the tip of a growing crack.

With the change of variables (5) and changing the variable of integration from $`\tau `$ to $`\zeta `$ gives

$$_{z_0}^{(\sqrt{z}\mathrm{})^2}𝑑\zeta \frac{(d\tau /d\zeta )}{\left[z\zeta \right]^{\frac{m}{2}}}+\frac{tt_{\mathrm{}}}{\left[(2\sqrt{z}\mathrm{})\mathrm{}\right]^{\frac{m}{2}}}=\left(\frac{3}{2\sigma _0}\right)^m\frac{d^{}}{z^{\frac{m}{2}}}\mathrm{with}\sqrt{z}a(t_{\mathrm{}})=\mathrm{}.$$
(51)
Note that $`tt_{\mathrm{}}`$ can also be written

$$\tau (z)\tau ((\sqrt{z}\mathrm{})^2)2\mathrm{}\sqrt{z}\frac{d\tau }{dz}|_z+𝒪\left(\mathrm{}^2\right).$$
(52)
The integral in the l.h.s. of (51) is analyzed similarly to the previous case (17). We multiply the integral by $`1/(yz)^{1\frac{m}{2}}`$ and integrate over $`z`$ from $`(\sqrt{z_0}+\mathrm{})^2`$ to $`y`$:

$$_{(\sqrt{z_0}+\mathrm{})^2}^y𝑑z\left(_{z_0}^{(\sqrt{z}\mathrm{})^2}𝑑\zeta \frac{(d\tau /d\zeta )}{\left(z\zeta \right)^{\frac{m}{2}}(yz)^{1\frac{m}{2}}}\right)=_{z_0}^{(\sqrt{y}\mathrm{})^2}𝑑\zeta \frac{d\tau }{d\zeta }\left[_{(\sqrt{\zeta }+\mathrm{})^2}^y\frac{dz}{\left(z\zeta \right)^{\frac{m}{2}}(yz)^{1\frac{m}{2}}}\right],$$
(53)
where we have used the equality of the triangle $`_{(\sqrt{z_0}+\mathrm{})^2}^y𝑑z_{z_0}^{(\sqrt{z}\mathrm{})^2}𝑑\zeta =_{z_0}^{(\sqrt{y}\mathrm{})^2}𝑑\zeta _{(\sqrt{\zeta }+\mathrm{})^2}^y𝑑z`$.

The integral in the bracket in the r.h.s. is the same as in (29), which defines the function $`b(m,ϵ)`$ with

$$ϵ(\zeta )=\frac{(\sqrt{\zeta }+\mathrm{})^2\zeta }{y\zeta }\frac{2\mathrm{}\sqrt{\zeta }}{y\zeta }.$$
(54)
Note that $`ϵ(\zeta )`$ is now a function of $`\zeta `$.

Using (52) and (53), expression (51) gives

$$_{z_0}^{(\sqrt{y}\mathrm{})^2}𝑑\zeta \frac{d\tau }{d\zeta }b(m,ϵ(\zeta ))+_{(\sqrt{z_0}+\mathrm{})^2}^y𝑑z\frac{2\mathrm{}\sqrt{z}}{(yz)^{1\frac{m}{2}}\left[(2\sqrt{z}\mathrm{})\mathrm{}\right]^{\frac{m}{2}}}\frac{d\tau }{dz}=\left(\frac{3\sqrt{d^{}}}{2\sigma _0}\right)^m_{(\sqrt{z_0}+\mathrm{})^2}^y\frac{d\zeta }{\zeta ^{\frac{m}{2}}(y\zeta )^{1\frac{m}{2}}}.$$
(55)
Since

$$b(m,ϵ(\zeta ))1/ϵ^{\frac{m}{2}1}\left(\frac{2\mathrm{}\sqrt{\zeta }}{y\zeta }\right)^{1\frac{m}{2}},$$
(56)
we see that the first integral of the l.h.s. of (55) is negligible compared to the second integral of the l.h.s. of (55) in the limit of large cracks, i.e. large $`z`$. Neglecting $`\mathrm{}`$ compared to $`\sqrt{z}`$ in the denominator of the integrant of the second integral of (55) and equating this second integral to the r.h.s. gives the following equation

$$\frac{d\tau }{dz}=(2\mathrm{})^{\frac{m}{2}1}\left(\frac{3\sqrt{d^{}}}{2\sigma _0}\right)^m\frac{1}{z^{\frac{m}{4}+\frac{1}{2}}}.$$
(57)
For $`m>2`$, $`\frac{m}{4}+\frac{1}{2}>1`$ and the solution of (57) is

$$a(t)\frac{\mathrm{}}{(t_ct)^{\frac{2}{m2}}},$$
(58)
where $`t_c`$ is determined from the initial size of the crack. The finite-time singularity results from the ever-increasing stress field at the fixed distance $`\mathrm{}`$ from the crack tip. This solution (58) is qualitatively different from the solution (24) found in the regime $`m<2`$ as (58) depends in a fundamental way upon the existence of the regularization scale $`\mathrm{}`$.

#### 5.2.2 Saturation by fixing an absolute maximum stress

An alternative prescription for the regularization is that the stress saturates at a constant value $`\sigma _{\mathrm{max}}`$. This is in constrast with the previous regularization scheme where the stress saturates at a value reached at a constant distance, this value thus increasing with the crack length. Expression (50) is then changed into

$$_0^t_{\mathrm{}}𝑑\tau \left(\frac{2\sigma _0}{3}\right)^m\frac{[a(t)]^m}{\left([a(t)]^2[a(\tau )]^2\right)^{\frac{m}{2}}}+[\sigma _{\mathrm{max}}]^m(tt_{\mathrm{}})=d^{},$$
(59)
with

$$\frac{2\sigma _0}{3}\frac{a(t)}{\left([a(t)]^2[a(t_{\mathrm{}})]^2\right)^{\frac{1}{2}}}=\sigma _{\mathrm{max}},$$
(60)
which is the condition that the stress saturates. It gives

$$[a(t_{\mathrm{}})]^2=[a(t)]^2(1A^2),$$
(61)
where

$$A=\frac{2\sigma _0}{3\sigma _{\mathrm{max}}}.$$
(62)
The equation (59) governing the dynamics of the crack tip can thus be written

$$\frac{2\sigma _0}{3}_{z_0}^{(1A^2)z}𝑑\zeta \frac{d\tau /d\zeta }{(z\zeta )^{m/2}}+[\sigma _{\mathrm{max}}]^mA^2z\frac{d\tau }{dz}=d^{},$$
(63)
where we have used the expansion

$$\tau (z)\tau ((1A^2)z)A^2z\frac{d\tau }{dz}|_z+𝒪\left(A^4\right),$$
(64)
valid in the interesting regime $`\sigma _{\mathrm{max}}>>\sigma _0`$ giving $`A<<1`$.

The second term in the l.h.s. of (63) dominates the first integral for large $`z`$, as can be checked a posteriori. For large crack sizes, the expression (64) can thus be simplified into

$$[\sigma _{\mathrm{max}}]^mA^2z\frac{d\tau }{dz}=d^{},$$
(65)
whose solution is

$$a(t)=a_0e^{t/t_0},$$
(66)
where

$$t_0=\frac{2d^{}}{A^2[\sigma _{\mathrm{max}}]^m}=\frac{9d^{}}{2\sigma _0^2[\sigma _{\mathrm{max}}]^{m2}}.$$
(67)
## 6 Beyond the mean field version by functional renormalization

Let us restrict our discussion to the case $`m=1`$ for which we have the complete analytical solution for the crack dynamics. The solution (15) with its asymptotic behavior (14) is not physically reasonable, as the crack reaches an infinite length in a finite time. The $`(t_ct)^1`$ singularity has been found to appear as the consequence of a geometric nonlinearity on an otherwise linearized mechanical problem. In reality, nonlinearity, viscosity, feedback, spatial heterogeneity of material properties and of cracking should modify the singularity. In addition, the main simplification in the previous approach is to neglect the impact of damage on the elastic coefficients of the material, thus leading to a stress field created by the crack which is identical to the field that the same static crack would generate in an undamaged material. Our hypothesis is that such modification can be deduced by a smooth or regular deformation of the solution previously obtained.

In this goal, we propose to apply the Yukalov-Gluzman functional renormalization method  to the series expansion of the solution (15) to obtain the renormalized law that accounts for these effects in a generic sense. Let us first consider the asymptotic power law singularity (14)

$$a(t)=\frac{2a_0}{\pi }(1x)^1,\mathrm{where}xt/t_c.$$
(68)
The powers $`x^n`$ in the expansion

$$a(t)=\frac{2a_0}{\pi }\left(1+x+x^2+x^3+\mathrm{}\right)$$
(69)
may be considered as hidden free parameters. Indeed, let us multiply the expansion by $`x^s`$. We then have a trial expansion for the solution. For $`s=0`$, we return to the regular expansion. Such multiplication can be applied repeatedly, for instance using the functional renormalization method . The idea behind the introduction of the multiplicative (control) function such as the power $`s`$ in $`x^s`$ is to deform smoothly the initial functional space of the expression $`a(t)`$ taken as an approximation to be improved. The condition for the improvement is to obtain a faster and better controlled convergence in the space of the modified functions upon addition of successive terms $`x^n`$ in the expansion. By this procedure, the dominant poles are eliminated or weakened as a result of a sequential reduction of stress level at each step of the resummation procedure. This corresponds to utilizing the information from the initial series pertaining to the times preceding the critical time $`t_c`$, where the level of damage is lower. Thus, the renormalization procedure is performing a mapping from the dynamics at early time far from the critical point to later times closer to the critical time. The stabilization stems from the fact that the information contained in the initial series related to times close to $`t_c`$ is minimized on the basis that it has an overly destabilizing effect in the description and should be weighted less than the information at earlier times.

At each step of the functional renormalization corresponding to the addition of a new term, we select the renormalized function according to the principle of minimum “local” multiplier, i.e., maximum stability on each sub-step of the renormalization procedure. Since these multipliers are proportional to the derivative $`da/dt`$ , the principle of minimal multiplier implies a selection of the real-time trajectory of the crack with minimal rate of damage (minimal stresses). In other words, this procedure amounts to improve the theory by allowing the crack to organize and develop so as to choose the most favorable path or dynamics. It can be shown  that, at each step, the choice of a formally infinite exponent $`s`$ corresponds to the minimal multiplier at arbitrary time.

The functional form of a super-exponential solution is selected by this procedure: starting from an expansion $`1+a_1x+a_2x^2+a_3x^3+\mathrm{}+a_kx^k`$, the renormalized expression is as follows. With the use of the notation

$$b_0=a_0,b_k=\frac{a_k}{a_{k1}},k=1,2,\mathrm{},$$
(70)
we obtain the bootstrap self-similar approximant up to order $`k`$

$$F_k(x)=b_0\mathrm{exp}(b_1x\mathrm{exp}(b_2x\mathrm{exp}(\mathrm{}b_{k1}x\mathrm{exp}(b_kx)))\mathrm{}),$$
(71)
introduced by Yukalov and Gluzman .

Let us now apply this result to the case (69) where all the coefficients $`a_n`$ are equal to $`1`$. The corresponding renormalized approximant replacing the initial input $`1/(1x)`$ of the expansion reads

$$F(x)=\mathrm{exp}\left(x\mathrm{exp}\left(x\mathrm{}\mathrm{exp}\left(x\right)\right)\mathrm{}\right).$$
(72)
This embedded exponential series converges to a well-defined function. To determine it, we note that $`F(x)`$ obeys the recursion relation

$$F_{k+1}(x)=\mathrm{exp}\left(xF_k(x)\right).$$
(73)
The fixed point to which these series of approximants converge is thus solution of

$$F=\mathrm{exp}\left(xF\right),$$
(74)
noting $`F=\pi a/2a_0`$. The limit $`F_{\mathrm{}}(x)`$ exists for $`\mathrm{e}x\mathrm{e}`$.

The fixed point $`F(x)`$ can be shown  to be the solution of the equation

$$\frac{dF}{dx}=\frac{F^2}{1xF}.$$
(75)
Searching for a solution in the form of a Taylor series

$$F(x)=\underset{n=0}{\overset{\mathrm{}}{}}y_nx^n,y_0=1,$$
(76)
we get

$$y_n=\frac{(n+1)^{n1}}{n!}.$$
(77)
Since $`n!n^ne^n`$, $`y_ne^n`$ for large $`n`$ and the generic term $`y_nx^n`$ in the series (76) is proportional to $`(ex)^n`$. This shows that the radius of convergence of the series (76) is $`1/e`$.

$`F(x)`$ has a singularity when $`x`$ approaches $`1/e`$ from below, whose shape is obtained by expansions of expression (74):

$$F(x)=_{x1/e}e\left(1\sqrt{2}e^{3/2}\sqrt{1/ex}\right).$$
(78)
Thus, the self-similar functional renormalization has transformed a pole (divergence of $`a`$) at $`t=t_c`$ into a square root singularity (finite $`a`$) at a smaller $`t=t_c/e`$. In this renormalized theory, the crack accelerates up to the time $`t_c/e`$ as which time its velocity diverges, while the crack is still finite. This announces the global breakdown. It is interesting that the exponent $`1/2`$ is close to the value found for acoustic emissions in experiments .

We can offer the following physical intuition for this transformation from the solution (14) with $`\beta =1`$ to to $`\beta =1/2`$. As the material becomes more and more damaged, the ulterior functional dependence of damage as a function of applied stress is modified. Actually, the series of functional renormalization amounts to effectively evolve or renormalize the damage law (2) into a succession of effective laws captured by the sequence of approximants, each approximant order corresponding to an increase in the overall damage of the material. Here, we have a mapping between a measure of evolution via the cumulative damage, i.e., a measure of passed time, and the order of the approximants and thus the distance to the fixed point in the functional space.

Consider now the general case $`a(t)(1x)^\beta `$. Expanding in power series, we get

$$(1x)^\beta =\underset{n=0}{\overset{\mathrm{}}{}}a_nx^n,$$
(79)
where

$$a_n=\frac{(n+\beta 1)!}{n!(\beta 1)!}.$$
(80)
The Yukalov-Gluzman renormalization scheme gives the superexponential (71) with coefficients $`b_n`$ given by

$$b_n=(n+\beta 1)/n.$$
(81)
Since $`b_n1`$ for large $`n`$ for any $`\beta `$, the fixed point of the approximants is controlled by the same finite square-root singularity of the type (78). Thus, the functional renormalization maps all finite-time singularities with different exponent $`\beta `$ on the same universal law $`a(t)=a(t_c)C\sqrt{t_ct}`$, where $`C`$ is a constant depending in particular on $`\beta `$.

## 7 Concluding remarks

Two main regimes have been found for the growth of a crack in a medium obeying the damage law $`d(d)/dt=\sigma ^m`$ (equation (2)), where $`\sigma `$ is the local stress. For $`0<m<2`$, a pre-existing crack grows to infinity in finite time and the divergence occurs as a power law finite-time singularity. For $`m2`$, the solution exists for all times but the characteristic time scale of the crack growth is an increasing function of a microscopic length scale, which is essential for regularizing the otherwise ill-defined problem. This microscopic length scale embodies the physical mechanism(s) by which the mathematical stress singularity at the crack tip of a perfectly sharp crack is rounded-off. We have examined two main scenarios, a damage-limited rupture and a stress-limited rupture.

The remarkable behavior of this simple model results from the form of the irreversible damage law, in particular from the fact that any non-vanishing stress increases the damage. Damage at any point is thus a kind of ever increasing counter of the history of the stress on that point. This feature prevents the existence of stationary solutions of cracks propagating at constant velocities. In contrast, we only obtain “run-aways.”

Stationary solutions can be obtained in simple generalizations of the damage law (2), for instance with a stress threshold below which no damage occurs or with a healing or work-hardening term allowing recovering of the material and decrease of the damage when the stress is low. Such situations have been investigated in discrete two-dimensional models .


format: text

text:
# A variational Monte Carlo calculation of dynamic multipole polarizabilities and van der Waals coefficients of the PsH system

## PACS number(s): 36.10.-k, 02.70.Lq

While experimentalists relay everyday on positrons and positronium atoms (Ps) to collect information about microscopic features of macroscopic systems like solutions, polymers and crystals, much less effort has been devoted to the theoretical understanding of the complex interactions that take place between ordinary matter and positrons. Among the explored avenues of this field, we mention the interest in predicting the stability of classes of compounds like e+M and MPs \[1-9\], where M represents an atomic or molecular system, and the calculation of the cross sections in the scattering process of e+ and Ps on a molecule or an atom \[10-16\].

On the contrary, the evaluation of the interaction energy between e+M or MPs and a molecule or atom is an almost unexplored issue . We believe this fact is primarily due to the need of a very accurate trial wave function to describe correctly the correlated motions of electrons and positrons. So far, only variational calculations with explicitly correlated Gaussians (ECG)  or Hylleraas-type functions \[7-9\], and the DMC method \[2-4\] have shown to be able to adequately recover the correlation energy in positron-containing systems.

Related to the calculation of the interaction energies is the calculation of second order properties of positron-containing systems, a problem whose surface has been barely scratched in the past . These properties, specifically the dynamic polarizabilities, are strictly related to the van der Waals coefficients that describe the long range interaction between systems , representing a way to tackle the problem of the asymptotic intermolecular interactions. Recently, Caffarel and Hess showed that these properties can be computed by means of quantum Monte Carlo simulations  connecting the imaginary-time-dependent dynamics of the unperturbed system with the transition probabilities of a reference diffusion process. In this work we apply a modified version of their method to compute dynamic multipole polarizabilities for PsH, H, He, and H- as a way to understand the behaviour of these systems when interacting with an external field, and as a first step towards the definition of the interaction potential between PsH and the ordinary matter.

As far as we know, the work by Le Sech and Silvi  is the only one reporting calculations on the effect of a constant electric field on PsH. In that work they computed both the static dipole polarizability, 123 a.u., and the behaviour of the annihilation rate $`\mathrm{\Gamma }_{2\gamma }`$ versus the intensity of the field employing explicitly correlated wave functions, numerical integration, and a variation-perturbation approach. As by-product of our calculations of the potential energy curve of the e+LiH system , we obtained an estimation of the static dipole polarizability of 49(2) a.u., a value quite different from the one computed by Le Sech and Silvi. Since we believe this difference to be too large to admit an explanation based on the different accuracy of the methods used to compute this value, we plan to solve this puzzle in this work.

In the method by Caffarel and Hess  the frequency-dependent second order correction to the ground state energy is written as a sum of the two time-centered autocorrelation functions of the perturbing potential $`V`$

$$E_\pm ^{(2)}(\omega )=_0^{\mathrm{}}e^{\pm t\omega }C_{VV}(t)𝑑t$$
(1)
where the autocorrelation function $`C_{VV}(t)`$ is given by

$$C_{VV}(t)=V(0)V(t)_{\mathrm{\Psi }_0^2}V(0)^2_{\mathrm{\Psi }_0^2}$$
(2)
Here, $`\mathrm{}_{\mathrm{\Psi }_0^2}`$ indicates that the average has to be taken using the Langevin dynamics that samples the square of the exact ground state wave function of the unperturbed system.

Caffarel and Hess  showed that it is possible to compute $`C_{VV}(t)`$ employing an optimized trial wave function and the pure-diffusion Monte Carlo (PDMC) method, an alternative algorithm to the commonly used DMC with branching, where each walker explicitly carries its own weight along all the simulation .

In their work on He and H2, Caffarel, Rérat, and Pouchan  reported that the autocorrelation function $`C_{VV}(t)`$ becomes dominated by the noise at large times, and this fact might be due to the fluctuations of the walker weights that increase during a PDMC simulation, while the value of the autocorrelation function itself becomes smaller. While the second effect is intrinsic to the stochastic method, the first can be reduced employing a more accurate trial wave function that is able to reduce the weight fluctuations. Another possibility, giving up the exactness of the method (i.e. not sampling the exact $`\mathrm{\Psi }_0^2`$), is represented by the sampling of a quite accurate trial wave function without carrying around the weight for each walker, a method we call Perturbation Theory variational Monte Carlo (PT-VMC). This algorithm can be useful for those systems whose autocorrelation function has a large decaying time, as the case of H- and PsH. This large decaying time will increase the fluctuations of the carried weights, and hence the statistical noise in the autocorrelation functions in the long-time region.

As a test of the correctness of our computer program and of the accuracy of the method, we computed the first three autocorrelation functions, and hence the dynamic polarizabilities up to the octupolar one, for the two systems H and He. The analytical forms of the perturbing potentials were taken from Ref. . While for H we employed the exact ground state wave function and compared with the analytic values of the multipole polarizability , for the He case we used a 25 term Hylleraas-type wave function optimized by means of the standard energy minimization . We fitted the numerical $`C_{VV}(t)`$ results of our simulations with a linear combination of three exponential functions

$$C_{VV}(t)\underset{i=1}{\overset{3}{}}a_ie^{\lambda _it}$$
(3)
in order to have an analytical representation of the autocorrelation functions at all the times. Since it is important to reproduce accurately the long time behaviour of $`C_{VV}`$, the smallest $`\lambda _i`$ in Eq. 3 was independently calculated fitting $`ln\left[C_{VV}\right]`$ in the long time region with a first order polynomial. This choice was found to improve sensibly the goodness of the total fitting in this time range.

These analytical representations of $`C_{VV}`$ allow us to compute easily the integrals in Eq. 1 and to obtain simple expressions of $`\alpha (\omega )`$. The parameters obtained by the fitting procedure are available from the authors upon requests.

For both systems we found excellent agreement of the static polarizabilities (H $`\alpha _{dip}=`$4.495 au, $`\alpha _{quad}=`$15.034 au, $`\alpha _{oct}=`$133.105 au; He $`\alpha _{dip}=`$1.382 au, $`\alpha _{quad}=`$2.401 au, $`\alpha _{oct}=`$10.367 au) with the exact results for H , with PDMC results by Caffarel, Rérat, and Pouchan , with Glover and Weinhold upper and lower bounds for He , and with the accurate results by Yan et al. .

At this point we would like to stress that, although in the PT-VMC method the walkers carry always a unitary weight because the branching process is absent, similarly to the PDMC method the time step has to be chosen short enough to produce only a small time step bias. For these two systems we found the time step of 0.01 hartree-1 to be adequate to compute statistically exact results.

As a check of the ability of the PT-VMC method to compute polarizabilities also for highly polarizable systems whose exact wave function is more diffuse than the one of He and H, we selected the hydride ion as test case. For this system we optimized a 5 term Hylleraas-type wave function whose average properties are shown in Table 1 together with the accurate results obtained in Ref. . Table 1 contains also the multipole static polarizabilities computed in this work employing a time step of 0.01 hartree-1, and the static polarizability computed by Glover and Weinhold . Comparing the mean values in Table 1, one can notice that our 5 term wave function gives lower values than the ones obtained in Ref.  except for $`r_{}`$. This fact may explain the underestimation of the $`\alpha _{dip}`$ by PT-VMC, that recovers 92(2) percent of the accurate value. Nevertheless, this result represents a fairly good estimation of the static dipole polarizability for H-, a quantity that appears difficult to compute even with more complex approaches .

As far as PsH is concerned, we computed the autocorrelation functions using two different trial wave functions, including one ($`\mathrm{\Psi }_T^1`$) and 28 ($`\mathrm{\Psi }_T^{28}`$) terms . The choice of two trial wave functions to guide the Langevin dynamics was aimed at testing the dependency of $`C_{VV}(t)`$ on the quality of the wave function itself.

Employing the PT-VMC method and our wave functions for PsH, we computed the autocorrelation functions for three perturbation potentials:

$`V_1=x_1+x_2x_p`$ (4)
$`V_2={\displaystyle \frac{3(x_1^2+x_2^2x_p^2)(r_1^2+r_2^2r_p^2)}{2}}`$ (5)
$`V_3=x_1^3+x_2^3x_p^3{\displaystyle \frac{3[x_1(y_1^2+z_1^2)+x_2(y_2^2+z_2^2)x_p(y_p^2+z_p^2)]}{2}}`$ (6)
where the subscripts 1 and 2 indicate the two electrons, while the subscript $`p`$ indicates the positron. These potentials are the cartesian forms of the dipole, quadrupole, and octupole moment operators for the PsH system. Figure 1 shows the averaged correlation functions for $`V_1`$, $`V_2`$, and $`V_3`$ as obtained by the VMC method employing the 28 term trial wave function. Each value of the correlation functions was computed employing roughly 1010 configurations. From Figure 1 one can note the effect at large evolution times of the dispersion of the “trajectories” used to compute the autocorrelation function. This effect makes difficult the reproduction of the long-time regime of these functions due to the exponential decay and the roughly constant statistical error introduced by the method. Moreover, the statistical error strongly depends on the perturbation potential whose autocorrelation function is computed, i.e., more specifically on the dispersion of its mean value over the $`\mathrm{\Psi }_T^2`$ distribution.

The results for the static multipole polarizabilities, i.e. for $`\omega =0`$, computed with both trial wave functions, are shown in Table 2. While for the dipole polarizabilities there is a good agreement between the two values, larger differences are present for the higher multipole polarizabilities. This fact is an indication of the different accuracy of the two functions in approximating the exact wave function at large distances from the nucleus. In fact it can be shown that if one approximates the autocorrelation functions taking care only of the excitation to the first state of the appropriate symmetry, the autocorrelation function is proportional to $`V_i^2V_i^2`$, where $`V_i`$ is the perturbing potential. Comparing the dipole results with the value obtained by Le Sech and Silvi , again the large difference of the computed polarizabilities is apparent. As a final check for this problem, we computed the energy of the PsH when immersed in a weak static electric field $`F`$ by means of standard DMC simulations adding the linear potential $`F(x_1+x_2x_p)`$. To make our simulations stable, i.e. to avoid the dissociation of the PsH, we truncated the effect of the linear potential at $`\left|x_i\right|=15`$ bohr. We fitted the DMC results by means of the simple polynomial $`a\alpha _{dip}F^2/2`$, where $`\alpha _{dip}`$ is the static dipole polarizability, obtaining $`\alpha _{dip}=`$ 42.3(8) a.u. We believe that this result, statistically indistinguishable from the $`\alpha _{dip}`$ obtained by the PT-VMC method, gives the definitive answer to the problem of the PsH polarizability. Nevertheless, the discrepancy between our PT-VMC and DMC $`\alpha _{dip}`$ and the one computed by Le Sech and Silvi  remains puzzling. In our experience , to compute the matrix elements they needed, millions of configurations must be used even for systems like PsH to avoid to be fooled by a false convergence. Unfortunately, Le Sech and Silvi did not report any information about the number of configurations they used to compute the integrals, so we cannot judge the numerical accuracy of their results.

An attempt to estimate the total accuracy of our $`\alpha `$ results can be made comparing the polarizability values obtained by the two wave functions. These differ by 10 percent at most, a value that we feel might give a conservative estimate of the relative errors for the higher multipolar fields.

As stated previously, although dynamical polarizabilities are interesting on their own, they represent the basis to compute van der Waals dispersion coefficients for the interaction between different systems. Therefore, following Ref. , we present the calculation of the $`C_6`$, $`C_8`$, and $`C_{10}`$ dispersion coefficients between H, He, and PsH as a first effort to obtain accurate information on the interaction between positronic systems and ordinary matter in the framework of the Born-Oppenheimer approximation and second-order perturbation theory.

Using the fitted parameters for H, He, and PsH we computed the coefficients for the interaction between the ordinary systems and between these and PsH. The values are reported in Table 3. Since the values for the H-H, H-He, and He-He coefficients are accurately known , we use them as a test of the accuracy of our approach: all the values differ from the accurate results by Yan et al.  at most by one part over hundreds.

Comparing the $`C_n`$’s for the ordinary systems with the ones for the interaction with PsH, it strikes that these last are more than an order of magnitude larger than the formers. These features, due to the larger PsH polarizability, indicates that positronic systems strongly interact with ordinary matter even at large distances. Unfortunately, nothing can be said about location and depth of the total potential minimum. This strongly depends also on the effect of the repulsion between the positron cloud and the H and He nuclei, so that we believe a supermolecule approach is needed. In a previous work  we computed the interaction energy between H and PsH, showing that this system could have a metastable state. Although the dispersion coefficients for the interaction between He and PsH are smaller than those for PsH and H, they might be large enough to give rise to a potential well that could support at least a stable state. If this turns out to be the case, the He-PsH system could be the lightest van der Waals (i.e. bound by means of dispersion forces) stable dimer.

ACKNOWLEDGMENTS

Financial support by the Universita’ degli Studi di Milano is gratefully acknowledged. The authors are indebted to the Centro CNR per lo Studio delle Relazioni tra Struttura e Reattivita’ Chimica for grants of computer time.

Figure captions:

Figure 1. Logarithm of the correlation functions of the perturbing potentials.

a This work (5 term wave function)
b Ref.
c Ref.

a Ref.


text:
# Statistical physics of adaptive correlation of agents in a market

## Acknowledgments

We are grateful to Andrea Cavagna, Irene Giardina and Matteo Marsili for useful comments and discussions. We acknowledge financial support from EPSRC Grant No. GR/M04426 and EC Grant No. ARG/B7-3011/94/27.


text:
# Nd elastic scattering as a tool to probe properties of 3N forces

## I Introduction

One major goal in nuclear physics is to establish the properties of nuclear forces and to understand nuclear phenomena by solving the many- nucleon Schrödinger equation driven by those elementary nuclear forces. Meson theory had an important impact for the construction of nuclear forces (both NN and 3N forces), but it lacks systematics like it would be given by an expansion parameter and the meson-nucleon vertices had to be parametrised in an ad-hoc manner. Nevertheless the one-pion exchange is undisputed and the most advanced formulation of meson-exchanges in the so called full Bonn potential  is remarkably successful. Because of its energy dependence - a consequence of deriving it by old fashioned time ordered perturbation theory - it is not useful in $`A>2`$ systems. Energy independent one-boson exchange versions thereof, however, are useful . In addition more phenomenological NN potentials have been constructed with the aim to describe the rich set of experimental NN data as precisely as possible. This leads to an often called new generation of realistic NN potentials: AV18 , CD Bonn , Nijm I, II and 93 . They describe the NN data set with an unprecedented precision of $`\chi ^2`$ per data point very close to one. Very recently an updated CD Bonn  appeared, which takes newest data into account but has not been used in the present article. An upcoming approach to construct nuclear forces in a systematic manner is chiral perturbation theory . First applications to three- and four-nucleon systems have been done .

In recent years it became possible to solve exactly three- and four-nucleon bound states using standard integration and differentiation methods . Stochastic techniques allow to go beyond A=4 and right now low energy states of nuclei up to A=8 are under control . In all cases those realistic NN forces fail to provide the experimental binding energies; there is clear cut underbinding. For instance this amounts to 0.5-1 MeV in case of three nucleons and to 2-4 MeV in case of 4He. A natural further step is the consideration of 3N forces. This is an even harder theoretical challenge and presently the most often used dynamical process is the $`\pi \pi `$ exchange between three nucleons with an intermediate excited nucleon state, the $`\mathrm{\Delta }`$ . This is augmented by further ingredients of various types as will be detailed below. By properly adjusting parameters one can achieve correct 3N and 4N binding energies and reaches even a fairly successful description of low energy bound states energies of up to A= 8 . However, in the latter case the results point to an insufficient spin-orbit splitting of nuclear levels in light nuclei as e.g. in 5He . This may be caused by a wrong spin structure of present day 3NF’s or by not well enough established $`{}_{}{}^{3}P_{j}^{}`$ NN force components. It will be interesting to see in the future the predictions of nuclear forces based on chiral perturbation theory.

Though this first signal on 3N force effects resulting from discrete states is important, a more detailed investigation of 3NF properties can be carried through in scattering processes, where a rich set of spin observables is available. The tremendous advance in computational resources allowed in recent years to make exact predictions for three- nucleon scattering using nuclear forces in all their complexities . Also experimentally one can access nowadays spin observables in Nd scattering where in the initial states the deuteron and / or the nucleon is polarized and after the reaction also the polarization of the outgoing particles can be measured . This leads to a very rich spectrum of observables in Nd elastic scattering and the Nd break-up processes. Such a set of spin observables will be a solid basis to test the 3N Hamiltonian. Using available model Hamiltonians one can provide guidance in selecting specific observables and energies which are most appropriate to see 3NF properties. It is the aim of this article to do exactly that and to compare the theoretical predictions with already existing data.

In section II we review briefly our 3N scattering formalism and display the 3NF model forces which are presently en vogue and which we use. Some technical details referring to the partial wave decomposition of the momentum space representation of the Urbana IX 3NF  are given in the Appendix. In this article we restrict ourselves to elastic Nd scattering and refer to a forthcoming article for the break-up process. Our predictions for various nuclear force combinations and the comparison to available data are given in section III. We conclude in section IV.

## II Scattering Formalism and 3NF Models

We refer to  for a general overview on 3N scattering and specifically our way to formulate it. For the inclusion of 3NF s we found meanwhile a more efficient way . It is a direct generalization of what is being used for the 3N bound state . We define an amplitude $`T`$ via our central equation

$$T=tP\varphi +(1+tG_0)V_4^{(1)}(1+P)\varphi +tPG_0T+(1+tG_0)V_4^{(1)}(1+P)G_0T$$
(1)
The initial channel state $`\varphi `$ occurring in the driving terms is composed of a deuteron and a momentum eigenstate of the projectile nucleon. The NN $`t`$-operator is denoted by $`t`$, the free 3N propagator by $`G_0`$ and $`P`$ is the sum of a cyclical and an anticyclical permutation of three particles. The 3N force $`V_4`$ can always be decomposed into a sum of three parts

$$V_4=V_4^{(1)}+V_4^{(2)}+V_4^{(3)},$$
(2)
where $`V_4^{(i)}`$ is symmetrical under the exchange of the nucleons $`jk`$ with $`jik`$. As seen in Eq. (1) only one of the three parts occurs explicitely, the others via the permutations contained in $`P`$. The physical break-up amplitude is given via

$$U_0=(1+P)T.$$
(3)
The Faddeev-like integral equation (1) has the nice property that its iteration inserted into Eq. (3) yields immediately the multiple scattering series, which gives a transparent insight into the reaction mechanism. Here in this article we concentrate on elastic scattering, whose amplitude is given by

$$U=PG_0^1+PT+V_4^{(1)}(1+P)\varphi +V_4^{(1)}(1+P)G_0T$$
(4)
The first term is the well known single particle exchange diagram, then there are terms where either $`V_4`$ or the $`t`$’s interact once and then the remaining parts result from rescattering among the three particles. Again inserting the iteration of $`T`$ as given in Eq. (1) into Eq. (4) yields a transparent insight .

The definition of the various spin observables can be found in . They have the general form

$$S^\mu _fI=\frac{1}{6}\underset{\nu }{}S^\nu _i\mathrm{Tr}\left(MS^\nu M^{}S^\mu \right),$$
(5)
where $`I`$ is the elastic cross section summed over the spin states in the final state, $`M`$ is the physical elastic scattering amplitude related directly to $`U`$ and $`S^\mu `$ is a suitable set of 3N spin operators.

We shall encounter nucleon and deuteron vector analyzing powers $`A_y(N)`$ and $`A_y(d)`$ (i$`T_{11}`$), where in the initial state either the nucleon (N) or the deuteron (d) is vector polarized. Further, the deuteron can be tensor polarized in the initial state leading to the three tensor analyzing powers $`T_{2k}`$ (k = 0,1,2). Also both particles can be polarized in the initial state leading to very many spin correlation coefficients $`C_{\alpha ,\beta }`$, where $`\alpha `$ refers to the spin directions of the nucleon and beta to vector and tensor polarizations of the deuteron. Further information on the dynamics can be found in spin transfer coefficients $`K_\alpha ^\beta ^{}`$, where $`\alpha `$ describes either a polarized nucleon or a polarized deuteron in the initial state and $`\beta ^{}`$ similarly the polarization for a particle in the final state. Of course all those quantities depend on the scattering angle.

Our nuclear model interaction consists of one of the NN forces mentioned in the introduction and a 3NF. For the 3NF we use the $`2\pi `$-exchange Tucson Melbourne (TM) model, a modified version thereof and the Urbana IX force. The TM model  has been around for quite some time. It is based on a low momentum expansion of the $`\pi N`$ off (the- mass) shell scattering amplitude. It has the form

$`V_4^{(1)}={\displaystyle \frac{1}{(2\pi )^6}}{\displaystyle \frac{g_{\pi NN}^2}{4m_N^2}}{\displaystyle \frac{\stackrel{}{\sigma }_2\stackrel{}{Q}}{\stackrel{}{Q}^{\mathrm{\hspace{0.17em}2}}+m_\pi ^2}}{\displaystyle \frac{\stackrel{}{\sigma }_3\stackrel{}{Q^{}}}{\stackrel{}{Q^{}}^{\mathrm{\hspace{0.17em}2}}+m_\pi ^2}}H\left(\stackrel{}{Q}^{\mathrm{\hspace{0.17em}2}}\right)H\left(\stackrel{}{Q^{}}^{\mathrm{\hspace{0.17em}2}}\right)`$

$$\left\{\stackrel{}{\tau }_2\stackrel{}{\tau }_3\left(a+b\stackrel{}{Q}\stackrel{}{Q^{}}+c(\stackrel{}{Q}^{\mathrm{\hspace{0.17em}2}}+\stackrel{}{Q^{}}^{\mathrm{\hspace{0.17em}2}})\right)+di\stackrel{}{\tau }_3\times \stackrel{}{\tau }_2\stackrel{}{\tau }_1\stackrel{}{\sigma }_1\stackrel{}{Q}\times \stackrel{}{Q^{}}\right\}.$$
(6)
The elements of the underlying Feynman diagram are obvious: the two pion propagators depending on the pion momenta $`\stackrel{}{Q}`$ and $`\stackrel{}{Q^{}}`$, the two $`\pi NN`$ vertex amplitudes and most importantly the parametrisation of the $`\pi N`$ amplitude inside the curly bracket which is combined with the isospins $`\stackrel{}{\tau }_2`$ and $`\stackrel{}{\tau }_3`$ of the two accompanying nucleons. On top of all that there is a strong form factor parametrisation given by

$$H\left(\stackrel{}{Q}^{\mathrm{\hspace{0.17em}2}}\right)=\left(\frac{\mathrm{\Lambda }^2m_\pi ^2}{\mathrm{\Lambda }^2+\stackrel{}{Q}^{\mathrm{\hspace{0.17em}2}}}\right)^2$$
(7)
In what we denote by the TM 3NF we use the original parameters $`a=1.13/m_\pi `$, $`b=2.58/m_\pi ^3`$, $`c=1.0/m_\pi ^3`$, $`d=0.753/m_\pi ^3`$. They incorporate among others the physics resulting from an intermediate $`\mathrm{\Delta }`$ in a static approximation. The cut-off parameter $`\mathrm{\Lambda }`$ is used to adjust the $`{}_{}{}^{3}H`$ binding energy separately for different NN forces . For the convenience of the reader we show the $`\mathrm{\Lambda }`$–values in Table 1. Of course in a meson exchange picture additional processes should be added containing other meson exchanges like $`\pi \rho `$, $`\rho \rho `$; also different intermediate excited states might play a role. To some extent 3NF models with respect to those extensions have already been developed and applied . Further studies should be performed.

The parametrisation of the TM 3NF has been criticized somewhat, since it violates chiral symmetry . A form consistent with chiral symmetry (though not a complete one to that order in the appropriate power counting) is obtained by modifying the $`c`$-term so that the long-range part is absorbed into the $`a`$-term, leading to a new $`a^{}a2m_\pi ^2c=0.87/m_\pi `$  what essentially means a change of sign for $`a`$ and that the short range part is dropped. This form will be called TM’ later on. The corresponding $`\mathrm{\Lambda }`$-value when used with the CD Bonn potential is $`\mathrm{\Lambda }=4.593m_\pi `$.

The two-meson exchange 3NF has also been studied by Robilotta et al.  leading to the Brazilian 3NF. It is similar to the one of TM and also the results gained for low energy Nd elastic scattering observables  are similar to the ones for the TM 3NF. In this article we do not take that force into account. Instead we included the Urbana IX 3NF , which is heavily used in the Urbana-Argonne collaboration. It will be interesting to see its effects for 3N scattering observables. At very low energies it has been used in that context before by the Pisa group . That force is based on the old Fujita-Mijazawa ansatz  of an intermediate $`\mathrm{\Delta }`$ occurring in the two-pion exchange and augmented by a spin independent short range piece. It has the form

$`V_4^{(1)}=A_{2\pi }\left[\{X_{12},X_{13}\}\{\stackrel{}{\tau }_1\stackrel{}{\tau }_2,\stackrel{}{\tau }_1\stackrel{}{\tau }_3\}+{\displaystyle \frac{1}{4}}[X_{12},X_{13}][\stackrel{}{\tau }_1\stackrel{}{\tau }_2,\stackrel{}{\tau }_1\stackrel{}{\tau }_3]\right]`$

$$+U_0T_\pi ^2(r_{12})T_\pi ^2(r_{13}),$$
(8)
where

$$X_{ij}=Y_\pi (r_{ij})\stackrel{}{\sigma }_i\stackrel{}{\sigma }_j+T_\pi (r_{ij})S_{ij},$$
(9)
with

$$Y_\pi (r)=\frac{e^{m_\pi r}}{m_\pi r}\left(1e^{cr^2}\right),$$
(10)
and

$$T_\pi (r)=\left[1+\frac{3}{m_\pi r}+\frac{3}{(m_\pi r)^2}\right]\frac{e^{m_\pi r}}{m_\pi r}\left(1e^{cr^2}\right)^2,$$
(11)
and where

$$S_{ij}=\mathrm{\hspace{0.17em}3}\stackrel{}{\sigma }_i\widehat{r}_{ij}\stackrel{}{\sigma }_j\widehat{r}_{ij}\stackrel{}{\sigma }_i\stackrel{}{\sigma }_j$$
(12)
is the tensor operator.

Since we work in momentum space and in a partial wave expansion the form given in Eq. (8) has to be rewritten. We could follow the steps laid our before for the corresponding representation of the TM 3NF  and delegate all that to the Appendix.

Since there is no apparent consistency of the mostly phenomenological realistic NN forces and the 3NF models we test various combinations thereof. In all cases, however, we require that the particular choice for the 2N interaction and the 3NF should reproduce the experimental triton binding energy. Some of the 3N observables scale with the triton binding energy . The adjustment to the triton binding energy has the advantage that our investigation is not misled by these scaling effects.

With respect to the intermediate $`\mathrm{\Delta }`$ one should say that very likely the static approximation is not justified and the $`\mathrm{\Delta }`$ should be allowed to propagate like the nucleon. This has been pursued intensively for instance by the Hannover group  and recent work has been also devoted to the 3N continuum .

In view of all that it is quite clear that our present study is not at all complete but can at most provide some insight, what kind of effects specific 3NF models might generate. As we shall see effects of that sort are needed, since NN force only predictions often fail to describe the data, especially at the higher energies. These challenges call for a systematic approach and at least for the leading spin structures chiral perturbation theory might be a good candidate . This is left to a future investigation. Here we concentrate on the current models and show their strengths and failures.

## III Predictions of 3NF Effects and Comparison to Data

Since we would like to cover a wide range of incoming neutron energies from below the nd breakup threshold up to 200 MeV it is necessary to take a sufficient number of partial wave states into account in order to get converged solutions of the Faddeev equations. In all calculations presented in this paper we went up to the two-nucleon subsystem total angular momentum $`j_{max}=5`$. This corresponds to a maximal number of $`142`$ partial wave states (often called channels) in the 3N system. We checked that the convergence has been achieved by looking at the results obtained for $`j_{max}=6`$, which increases the number of channels to 194. This convergence check refers to a calculation without a 3NF. The inclusion of 3NF’s has been carried through for all total angular momenta of the 3N system up to $`J=13/2`$. These high angular momenta are required at the higher energies $`100`$ MeV. The longer ranged 2N interactions require states up to $`J=25/2`$ at the higher energies in order to get converged results.

A phenomenological criterium for 3NF effects is that the data lie outside the spread of NN force predictions only. In the following figures we shall always include a shaded band (called “band1”), which covers the predictions of the AV18, CD-Bonn, NijmI, II and 93 NN forces. Unfortunately we cannot include the pp-Coulomb force in our approach and thus have to live with some theoretical uncertainty when comparing to pd-data. At the higher energies, however, those effects should be small. Also in case of AV18 we do not take the various electromagnetic corrections into account, which leads for example to a slightly wrong deuteron binding energy ($`E_d`$ = 2.242 instead of 2.225 MeV). This, however, has only a small effect on our results, which is mostly of kinematical origin, since the phase-shifts obtained without those additional terms differ only slightly from the standard ones. The kinematical effects are seen predominantly in the breakup process, where this small defect in the deuteron binding energy leads to correspondingly small shifts in peak structures like for instance final state interaction peaks.

We shall combine another group of curves into a band. The TM 3NF can be combined with the five NN forces. In all cases the cut-off value $`\mathrm{\Lambda }`$ in Eq. (7) has been adjusted separately for each NN force to the 3H binding energy . Since that interplay is a purely phenomenological step the outcome is theoretically not under control and we combine all the results into a second band (called “band2” in the following). Next we want to compare the TM 3NF and the modified TM’, which is more consistent with chiral symmetry. We combine it with CD-Bonn and show the CD-Bonn+TM’ prediction as a dashed curve. Finally we compare the TM and the Urbana IX 3NF’s and combine them with AV18. The combination AV18+URBANA IX will appear as a solid curve. There are clearly more combinations possible but it is sufficient to get an orientation on the magnitudes of expected effects.

We begin with the differential cross section in Fig. 1. The various NN force predictions are rather close together with a small spread in the minima. Including the TM 3NF there is again a small spread in the minima (practically negligible at 3 MeV) but the minima are shifted upwards, rather well into the data  except at 3 MeV, where the 3NF prediction is shifted slightly downwards. The phenomenon of shrinkage of the spread between different NN potential predictions by including a 3NF is often called a scaling phenomenon. It occurs at low energies and is related to the three-nucleon binding energy, which by construction is common to all those curves in band2. The TM and TM’ 3NF’s together with CD-Bonn give slightly different predictions in the minima especially at the two highest energies. (The CD Bonn+TM prediction lying inside band2 is not shown.) In the backward angular region they differ significantly and the 135 MeV precise backward angular distributions data prefer the TM 3NF. (See insertion in Fig. 1; again the CD Bonn+TM prediction is not explicitely shown.) On the other hand TM and URBANA IX together with AV18 are very similar at the two higher energies but differ significantly at $`E`$= 65 MeV. Certainly Coulomb force effects should be taken into account at this energy before a final conclusion can be drawn. The few nd data near the minimum would strongly disagree with all our 3NF predictions and a confirmation (or rejection) would be highly desirable. However, independent from possible Coulomb force contributions, it is clearly seen, that even such a simple observable as the elastic scattering differential cross section exhibits large 3NF effects at higher energies. These effects are not trivial and depend not only on the incoming energy and the angle but also on the particular 3NF used. This calls for precise data for this observable in order to study the 3NF properties.

Let us now regard a selection out of the many spin observables in elastic Nd scattering. Fig.2 shows $`A_y(N)`$. The band for NN force predictions is always rather narrow, whereas band2 for the lowest and highest energy is distinctly broader. The two bands are separated predicting clearly 3NF effects especially at higher energies. The TM and TM’ predictions are distinctly different as well as the TM and Urbana IX predictions. It is interesting to note that here TM’ with CD Bonn and Urbana IX with AV18 are very similar (except at the lowest energy) and predict only small 3NF effects at 65 MeV, which are compatible with the $`A_y`$ data at this energy. At higher energies their effects become quite different from the TM ones. While in the region of the $`A_y`$ minimum around $`\theta _{cm}110^o`$ they increase $`A_y`$ as compared to the pure 2N force predictions, their action decreases $`A_y`$ in the backward angular region, contrary to the action of the TM 3NF model, bringing the theory closer to the data. In this way the TM’ and Urbana IX 3NF’s seem to solve partially the $`A_y`$ problem found at higher energies in . At 3 MeV the clear discrepancy of all theoretical predictions to the nd data of  is seen. At such a low energy it is well known that Coulomb force effects are large for $`A_y`$ decreasing significantly its maximum when compared to nd data . Thus also pd data lie very clearly (due to much smaller error bars) above all theoretical predictions. We see that this very well known low energy $`A_y`$ puzzle  cannot be solved by the 3NF models we are using in this article. A slightly increased maximum of $`A_y`$ for TM’ is far too small to play any significant role and possibly the solution should be also sought in an improvement of the $`{}_{}{}^{3}P_{j}^{}`$ NN force components  to which low energy $`A_y`$ is very sensitive. We would also like to point to the very recent result based on chiral perturbation theory , where $`A_y`$ can be described quite well in next-to-leading order (NLO). In that order of the power counting 3NF’s do not yet contribute. Those effective chiral NN forces are very different from the conventional ones. But this NLO result is just an intermediate step and the final answer has to wait for higher order contributions, which improve systematically the observables in the 2N, 3N, … systems at the same time.

For $`iT_{11}`$ shown in Fig. 3 the two bands are distinctly different and clearly the TM-band is supported by the data at the higher energies. In that case TM and TM’ are close together and also TM and URBANA IX except around 120  at the highest energy. The data shown in the figure for 190 MeV are taken at 197 MeV. Unfortunately they are absent around 120 .

Next we regard the three tensor analyzing powers $`T_{20}`$, $`T_{21}`$ and $`T_{22}`$ in Figs. 4-6. For $`T_{20}`$ the situation is very challenging. At 135 MeV the data between about 120 -150  do not agree with the overlapping bands and above 150  they lie just between the two bands. At small angles they agree with the overlapping bands and follow then the NN force prediction. In addition at the higher energies TM and TM’ differ as well as TM and URBANA IX. The strong deviations of the theory to the data at 3 MeV is simply caused by Coulomb force effects .

For $`T_{21}`$ the situation is different. The two bands are clearly distinct. Again the different 3NF’s predictions deviate strongly from each other. The data at 135 MeV follow more band1 than band2 and at the small angles TM’ or URBANA IX are preferred. Clearly this is a rather contradictory situation.

For $`T_{22}`$ the bands differ but the special 3NF predictions (dashed and solid lines) are similar but lie outside band2. Except at very backward and forward angles, where all curves essentially coincide, there is disagreement with the data at $`E`$= 135 MeV. For all three $`T_{2k}`$’s data at 190 MeV would be very valuable, since the various theoretical predictions differ dramatically.

There are many spin-transfer coefficients and we selected more or less arbitrarily five of them. In Figs. 7-11 we show $`K_{yz}^x^{}`$, $`K_{yy}^y^{}`$, $`K_{xx}^y^{}`$, $`K_{xz}^y^{}`$ and $`K_x^{y^{}z^{}}`$. For $`K_{yz}^x^{}`$ the bands strongly deviate at the higher energies and TM and TM’ as well as TM and URBANA IX differ drastically. The deviation of the bands from each other is less pronounced for $`K_{yy}^y^{}`$ and also the different 3NF predictions are less distinct except at the highest energy. Two data points agree with the 3NF predictions, while one, at 150 , is below any theoretical prediction. For $`K_{xx}^y^{}`$ the bands differ only at the lowest energy ($`E`$= 3 MeV) significantly. Otherwise the 3NF effects are rather modest and also differences among the different 3NF’s. The data at 135 MeV go across the various predictions. For $`K_{xz}^y^{}`$ the bands differ at the high energies and the two special 3NF predictions deviate significantly from the TM predictions. The one data point at 135 MeV lies somewhere in between. Finally $`K_x^{y^{}z^{}}`$ show again dramatic effects in relation to the two different bands. Also the special 3NF predictions are clearly different at the two higher energies.

Lastly we regard in Figs. 12-14 three different spin correlation coefficients $`C_{xy,x}`$, $`C_{yy}`$ and $`C_{zz}`$. The effects are dramatic for $`C_{xy,x}`$: the bands are quite different and also the special 3NF predictions. For $`C_{yy}`$ there are data  at $`E_p`$ = 197 MeV, which we inserted into the figure for 190 MeV. The data support the curves inside band2. The AV18+URBANA IX is significantly closer to experiment than CD Bonn+TM’. Finally for $`C_{zz}`$ the bands differ a lot at the high energies and the special 3NF predictions differ from the TM ones as well. One notes that URBANA IX and TM’ are similar.

## IV Summary and Conclusions

We performed a study of popular present day 3NF models with respect to the effects they cause in the 3N continuum. Based on the comparison of the realistic NN force predictions alone (“band1”) to the predictions of all NN forces combined with the TM 3NF (“band2”) one sees in many spin observables very drastic effects, which should clearly be discernible by experiments. On top the three different 3NF models, TM, TM’ and URBANA IX combined (arbitrarily) with CD-Bonn and AV18 lead again to other very distinct predictions. Specifically the effects are angular and energy dependent becoming large at higher energies. It seems that with a sufficiently rich and precise data basis such diversity of effects should allow to nail down the proper spin structure of 3NF’s.

Unfortunately there are up to now only few data available. The ones for the differential cross section support the shift in theory caused by 3NF’s. The existing high energy cross section data in the backward angular region prefer the structure of the TM 3NF. Also the existing deuteron vector analyzing powers at higher energies support rather well the predicted 3NF effects. On the other hand this three-body interaction predicts too large effects for the nucleon analyzing power. This observable seems to prefer the modified version of the Tucson-Melbourne model, TM’, which is consistent with chiral symmetry or the URBANA IX. For the tensor analyzing powers the situation is totally chaotic, for some scenarios one finds agreement, for others a strong disagreement: there is no preference for any of them. Clearly we are at the very beginning in investigating the spin-structure of the 3NF. Nevertheless the effects of all those 3NF’s are typically of the right order in magnitude, when they can be checked against data but the signs are not yet under control. The spin transfer coefficients carry also a lot of information and in some of them the two bands differ very much. Finally, the spin correlation coefficients appear also very informative and the very first data for $`C_{yy}`$ support the TM and URBANA IX.

Altogether the only conclusion possible is that the most popular current 3NF models show a lot of effects and data are needed to provide constraints. There is hope that further theoretical work guided by the chiral effective field theory approach will help to establish the proper spin structure of the three-nucleon force.

###### Acknowledgements.

This work was supported by the Deutsche Forschungsgemeinschaft (H.K., A.N. and J.G.), the Polish Committee for Scientific Research under Grant No. 2P03B02818 and the Science and Technology Cooperation Germany-Poland. W.G. would like to thank the Foundation for Polish Science for the financial support during his stay in Cracow. The numerical calculations have been performed on the Cray T90 and T3E of the NIC in Jülich, Germany.

## A Partial wave decomposition of the Urbana IX 3NF in momentum space.

The Urbana 3NF in Eq.(8) has to be put into a form suitable for the evaluation in partial wave decomposition. Therefore we rewrite Eq.(8) using the (anti)commutator of the isospin operators

$`\{\stackrel{}{\tau }_1\stackrel{}{\tau }_2,\stackrel{}{\tau }_1\stackrel{}{\tau }_3\}`$ $`=`$ $`2\stackrel{}{\tau }_2\stackrel{}{\tau }_3`$ (A1)
$`[\stackrel{}{\tau }_1\stackrel{}{\tau }_2,\stackrel{}{\tau }_1\stackrel{}{\tau }_3]`$ $`=`$ $`2i\stackrel{}{\tau }_1\stackrel{}{\tau }_2\times \stackrel{}{\tau }_3`$ (A2)
$`V_4^{(1)}`$ reads then

$`V_4^{(1)}`$ $`=`$ $`A_{2\pi }2\left(\stackrel{}{\tau _2}\stackrel{}{\tau _3}+{\displaystyle \frac{i}{4}}\stackrel{}{\tau _1}\stackrel{}{\tau _2}\times \stackrel{}{\tau _3}\right)X_{12}X_{13}`$ (A5)
$`+A_{2\pi }2\left(\stackrel{}{\tau _2}\stackrel{}{\tau _3}{\displaystyle \frac{i}{4}}\stackrel{}{\tau _1}\stackrel{}{\tau _2}\times \stackrel{}{\tau _3}\right)X_{13}X_{12}`$
$`+U_0T_\pi ^2(r_{12})T_\pi ^2(r_{13})`$
$`=`$ $`A_{2\pi }2\left(\stackrel{}{\tau _2}\stackrel{}{\tau _3}+{\displaystyle \frac{i}{4}}\stackrel{}{\tau _1}\stackrel{}{\tau _2}\times \stackrel{}{\tau _3}\right)X_{12}X_{13}`$ (A8)
$`+A_{2\pi }2\left(\stackrel{}{\tau _2}\stackrel{}{\tau _3}{\displaystyle \frac{i}{4}}\stackrel{}{\tau _1}\stackrel{}{\tau _2}\times \stackrel{}{\tau _3}\right)X_{13}X_{12}`$
$`+{\displaystyle \frac{U_0}{2}}\left(T_\pi ^2(r_{12})T_\pi ^2(r_{13})+T_\pi ^2(r_{13})T_\pi ^2(r_{12})\right)`$
Because the $`T^2`$ operators commute which each other, we could choose the symmetrical form in the last line.

The aim is to find matrix elements with respect to our standard basis states

$$|pq\alpha >_i|pq(ls)j(\lambda \frac{1}{2})J(jJ)𝒥(t\frac{1}{2})TM_T>_i$$
(A9)
Here the magnitudes of the Jacobi momenta of the subsystem and the outer particle are $`p`$ and $`q`$, respectively. The angular dependence is expanded in partial waves. Corresponding to $`p`$ and $`q`$ we introduce angular momenta $`l`$ and $`\lambda `$. As indicated in Eq.(A9) the orbital angular momenta couple with the spin of the subsystem $`s`$ and the outer particle $`\frac{1}{2}`$ to the total spin of the subsystem $`j`$ and outer particle $`J`$. These angular momenta are combined to the total spin $`𝒥`$ and its third component $``$. The total isospin of the subsystem $`t`$ couples with the isospin of the outer particle $`\frac{1}{2}`$ to $`T`$ and its third component $`M_T`$. Because we will make use of several sets of Jacobi momenta, we append an index $`i`$ which gives the number of the outer particle.

According to Eqs.(1) and (4) the operator $`V_4^{(1)}`$ acts on completely antisymmetric 3N states $`|\chi >=(1+P)|\varphi >`$ or $`|\chi >=(1+P)G_0T`$. In the next paragraphs we would like to establish some consequences of this fact.

To that aim we introduce a short hand notation for our basis states:

$$|(jk)i>|pq\alpha >_i$$
(A10)
In contrast to Eq.(A9) this definition fixes the ordering within the subsystem. $`(ij)`$ is understood to fix the momentum vector to $`\stackrel{}{p}=\frac{1}{2}(\stackrel{}{k}_i\stackrel{}{k}_j)`$ and the spin and isospin coupling within the subsystem to $`(s_is_j)s`$ and $`(t_it_j)t`$. In consequence we can distinguish $`|(ij)k>`$ and $`|(ji)k>`$ in this notation. In the partial wave decomposition there is a simple phase relation connecting both states

$$|(ij)k>=()^{l+s+t}|(ji)k>()^{(ij)}|(ji)k>$$
(A11)
Of course all sets of basis states are complete, therefore we can expand the incoming state in several ways

$`|\chi >`$ $`=`$ $`{\displaystyle |(31)2><(31)2|\chi >}`$ (A12)
$`=`$ $`{\displaystyle |(12)3><(12)3|\chi >}`$ (A13)
Due to the total antisymmetry of $`|\chi >`$, its matrix elements are equal

$$<(31)2|\chi >=<(12)3|\chi >=<(23)1|\chi >$$
(A14)
In consequence we can expand the different terms in Eq.(A5) on the right hand side in different Jacobi coordinates. Let us begin with the two $`U_0`$ terms in Eq.(A5):

$`M_3`$ $`=`$ $`{\displaystyle \frac{1}{2}}<(23)1|(12)3^{\prime \prime }><(12)3^{\prime \prime }|T_\pi ^2(r_{12})|(12)3^{\prime \prime \prime }>`$ (A18)
$`<(12)3^{\prime \prime \prime }|(31)2^{\prime \prime \prime \prime }><(31)2^{\prime \prime \prime \prime }|T_\pi ^2(r_{13})|(31)2^{}>`$
$`+{\displaystyle \frac{1}{2}}<(23)1|(31)2^{\prime \prime }><(31)2^{\prime \prime }|T_\pi ^2(r_{13})|(31)2^{\prime \prime \prime }>`$
$`<(31)2^{\prime \prime \prime }|(12)3^{\prime \prime \prime \prime }><(12)3^{\prime \prime \prime \prime }|T_\pi ^2(r_{12})|(12)3^{}>`$
In this equation we introduced several completeness relations and projected on two kinds of incoming states keeping in mind that one can add up the matrix elements because of the total antisymmetry of $`|\chi >`$. In the way we inserted the completeness relations in Eq.(A18), the matrix elements of $`T_\pi ^2`$ are evaluated in their natural coordinates. In this form the operator is diagonal in the quantum numbers and momenta of the outer particle and additionally it conserves the symmetry with respect to the interchange of the particles of the subsystem.

The matrix element $`<(12)3|T_\pi ^2(r_{12})|(12)3^{}>`$ depends on Jacobi coordinates and quantum numbers which single out the subsystem (12). By renumbering the particles one finds

$$<(12)3|T_\pi ^2(r_{12})|(12)3^{}>=<(31)2|T_\pi ^2(r_{13})|(31)2^{}>$$
(A19)
Note that the matrix element on the right hand side depends on coordinates which differ from the ones on the left hand side. Though the matrix elements are equal, one has to keep in mind that the actual meaning of the momenta and quantum numbers is different on both sides.

In the same manner one can find an important relation between the cyclic and anti cyclic transformations

$`<(ij)k|(jk)i^{}>=<(ik)j|(kj)i^{}>=()^{(ik)}()^{(kj)^{}}<(ki)j|(jk)i^{}>`$ (A20)
We would like to emphasize that the transformation itself doesn’t conserve the symmetry of the subsystem and therefore the completeness relation in ′′ and ′′′ states in Eq.(A18) have to include also the unphysical symmetric states in the two-body subsystems. With the help of Eq.(A20) and renumbering the particles in the second part of Eq.(A18) one finds

$`M_3`$ $`=`$ $`{\displaystyle \frac{1}{2}}<(23)1|(12)3^{\prime \prime }><(12)3^{\prime \prime }|T_\pi ^2(r_{12})|(12)3^{\prime \prime \prime }>`$ (A24)
$`<(12)3^{\prime \prime \prime }|(31)2^{\prime \prime \prime \prime }><(31)2^{\prime \prime \prime \prime }|T_\pi ^2(r_{13})|(31)2^{}>`$
$`+{\displaystyle \frac{1}{2}}()^{(23)}()^{(12)^{\prime \prime }}()^{(12)^{\prime \prime \prime }}()^{(31)^{\prime \prime \prime \prime }}<(23)1|(12)3^{\prime \prime }><(12)3^{\prime \prime }|T_\pi ^2(r_{12})|(12)3^{\prime \prime \prime }>`$
$`<(12)3^{\prime \prime \prime }|(31)2^{\prime \prime \prime \prime }><(31)2^{\prime \prime \prime \prime }|T_\pi ^2(r_{13})|(31)2^{}>`$
Because $`T_\pi ^2`$ conserves the symmetry of the subsystem $`()^{(12)^{\prime \prime }}=()^{(12)^{\prime \prime \prime }}`$ and $`()^{(31)^{\prime \prime \prime \prime }}=()^{(31)^{}}`$. Therefore Eq.(A24) reduces to

$`M_3`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(1+()^{(23)}()^{(31)^{}}\right)<(23)1|(12)3^{\prime \prime }><(12)3^{\prime \prime }|T_\pi ^2(r_{12})|(12)3^{\prime \prime \prime }>`$ (A26)
$`<(12)3^{\prime \prime \prime }|(31)2^{\prime \prime \prime \prime }><(31)2^{\prime \prime \prime \prime }|T_\pi ^2(r_{13})|(31)2^{}>`$
The incoming state is antisymmetric in the subsystem (31) hence $`()^{(31)^{}}=1`$. Therefore $`M_3`$ is zero for outgoing states which are symmetric in the subsystem (23) and equals twice the first part for the antisymmetric outgoing states. We restrict the outgoing states to antisymmetric ones. Then it is justified to write

$`M_3`$ $`=`$ $`<(23)1|(12)3^{\prime \prime }><(12)3^{\prime \prime }|T_\pi ^2(r_{12})|(12)3^{\prime \prime \prime }>`$ (A28)
$`<(12)3^{\prime \prime \prime }|(31)2^{\prime \prime \prime \prime }><(31)2^{\prime \prime \prime \prime }|T_\pi ^2(r_{13})|(31)2^{}>`$
We would like emphasize again that we restrict the outgoing, incoming and ′′′′ states to physical, antisymmetric states in the subsystems. Due to the coordinate transformations, the ′′ and ′′′ sums are not restricted anymore to antisymmetric subsystem states and the completeness relations run over all symmetries.

Lets turn to the $`A_{2\pi }`$ parts now. In the same manner we define the matrix elements of the $`A_{2\pi }`$ parts

$`\stackrel{~}{M}_3`$ $`=`$ $`<(23)1|(12)3^{\prime \prime }><(12)3^{\prime \prime }|X_{12}|(12)3^{\prime \prime \prime }>`$ (A32)
$`<(12)3^{\prime \prime \prime }|2\left(\stackrel{}{\tau _2}\stackrel{}{\tau _3}+{\displaystyle \frac{i}{4}}\stackrel{}{\tau _1}\stackrel{}{\tau _2}\times \stackrel{}{\tau _3}\right)|(31)2^{\prime \prime \prime \prime }><(31)2^{\prime \prime \prime \prime }|X_{13}|(31)2^{}>`$
$`+<(23)1|(31)2^{\prime \prime }><(31)2^{\prime \prime }|X_{13}|(31)2^{\prime \prime \prime }>`$
$`<(31)2^{\prime \prime \prime }|2\left(\stackrel{}{\tau _2}\stackrel{}{\tau _3}{\displaystyle \frac{i}{4}}\stackrel{}{\tau _1}\stackrel{}{\tau _2}\times \stackrel{}{\tau _3}\right)|(12)3^{\prime \prime \prime \prime }><(12)3^{\prime \prime \prime \prime }|X_{12}|(12)3^{}>`$
The isospin operators commute with the $`X_{ij}`$ operators. We combined them with the inner transformation matrix. Because $`\stackrel{}{\tau _2}\stackrel{}{\tau _3}`$ is symmetric with respect to the interchange of particle 2 and 3 and $`\stackrel{}{\tau _1}\stackrel{}{\tau _2}\times \stackrel{}{\tau _3}`$ is antisymmetric, the steps leading to Eq.(A20) can be repeated with the isospin matrix element. In addition there is a sign change for the antisymmetric part of the isospin operators.

$`<(12)3^{\prime \prime \prime }|2\left(\stackrel{}{\tau _2}\stackrel{}{\tau _3}+{\displaystyle \frac{i}{4}}\stackrel{}{\tau _1}\stackrel{}{\tau _2}\times \stackrel{}{\tau _3}\right)|(31)2^{\prime \prime \prime \prime }>`$ (A33)
$`=()^{(31)^{\prime \prime \prime }}()^{(12)^{\prime \prime \prime \prime }}<(31)2^{\prime \prime \prime }|2\left(\stackrel{}{\tau _2}\stackrel{}{\tau _3}{\displaystyle \frac{i}{4}}\stackrel{}{\tau _1}\stackrel{}{\tau _2}\times \stackrel{}{\tau _3}\right)|(12)3^{\prime \prime \prime \prime }>`$ (A34)
The symmetry of the $`X_{ij}`$’s leads to equal phases on both sides of their matrix elements: $`()^{(31)^{\prime \prime \prime \prime }}=()^{(31)^{}}`$ and $`()^{(12)^{\prime \prime }}=()^{(12)^{\prime \prime \prime }}`$. The matrix element of the $`A_{2\pi }`$ parts reduces to

$`\stackrel{~}{M}_3`$ $`=`$ $`\left(1+()^{(23)}()^{(31)^{}}\right)<(23)1|(12)3^{\prime \prime }><(12)3^{\prime \prime }|X_{12}|(12)3^{\prime \prime \prime }>`$ (A36)
$`<(12)3^{\prime \prime \prime }|2\left(\stackrel{}{\tau _2}\stackrel{}{\tau _3}+{\displaystyle \frac{i}{4}}\stackrel{}{\tau _1}\stackrel{}{\tau _2}\times \stackrel{}{\tau _3}\right)|(31)2^{\prime \prime \prime \prime }><(31)2^{\prime \prime \prime \prime }|X_{13}|(31)2^{}>`$
and its is again justified to restrict the outgoing states to antisymmetric subsystems and write

$`\stackrel{~}{M}_3`$ $`=`$ $`2<(23)1|(12)3^{\prime \prime }><(12)3^{\prime \prime }|X_{12}|(12)3^{\prime \prime \prime }>`$ (A38)
$`<(12)3^{\prime \prime \prime }|2\left(\stackrel{}{\tau _2}\stackrel{}{\tau _3}+{\displaystyle \frac{i}{4}}\stackrel{}{\tau _1}\stackrel{}{\tau _2}\times \stackrel{}{\tau _3}\right)|(31)2^{\prime \prime \prime \prime }><(31)2^{\prime \prime \prime \prime }|X_{13}|(31)2^{}>`$
It is standard to work out the two-body partial wave matrix elements for $`X_{ij}`$ and $`T_\pi ^2`$ in terms of Bessel transforms. Also the isospin and coordinate transformation matrix elements are standard and we refer to .