text:
# I Introduction

## I Introduction

Manipulating the outcome of quantum dynamics phenomena by a properly tailored external control is a topic of increasing activity \[1-15\]. Problems where the external control is electromagnetic have received the most attention, although other applications arise as well. Various implementations of quantum control have been experimentally realized \[16-25\]. A variety of control strategies have been suggested, and optimal control theory(OCT) provides the most general framework for acheiving field designs. Such designs will generally require further refinement in the laboratory through learning techniques to overcome design uncertainities \[17-20\]. A basic difficulty in attaining the control designs is the computational effort called for in solving the time-dependent Schroedinger equation, often repeatedly in an iterative fashion. This paper introduces the quantum fluid dynamic(QFD) control formulation to simplify this task. Calculations have shown  that QFD is capable of being much more efficient than conventional methods(e.g., FFT propagation), and this savings should carry over to the control design task. This paper will show how OCT can be combined with QFD.

The theoretical basis for dynamic control \[2-13\] is to create non-stationary states of one’s choice, by optimally designing the control field. Typically, the problem is posed as seeking an optimal field to drive a quantum wave packet to a desired target at a chosen time t=T. In the traditional approach \[3-10\] to quantum optimal control an objective design functional $`\overline{J}`$ is defined, which depends on the system wave function, a wave function like Lagrange multiplier, and the external field. Minimization of the objective functional leads to the identification of external field(s) capable of delivering a specific outcome. This process requires solving for the complex oscillatory wave function and the similarly behaved Lagrange multiplier, and due care is needed for their proper representation in a suitable basis for capturing their behaviour. Often the rapidly varying spatio-temporal behaviour of these functions necesitates the use of many unknowns in the basis.

This paper explores an alternative formulation for OCT to design the electric field. The formulation is based on the fluid dynamic view point of quantum mechanics \[26-29\], which by-passes the typically oscillatory nature of the wave function to exploit the generally smooth behaviour of the real density and the quantum current variables. Recent illustrations have demonstrated the smooth spatial and temporal nature of the variables and the ability to discretize them on a relatively small number of grid points . As background in section 2 we give a brief summary of the QFD formulation of quantum mechanics. Section 3 presents OCT within the framework of QFD for designing an electric field to meet a specific objective. Section 4 applies the OCT-QFD formulation for the manipulation of HCl. Section 5 concludes the paper.

## II Quantum Fluid Dynamics

The treatment below considers a single particle of reduced mass m, but the QFD formulation has an immediate extension to many particles. The time-dependent Schrödinger equation is given by

$`[{\displaystyle \frac{\mathrm{}^2}{2m}}^2+V+V_{ext}]\mathrm{\Psi }(𝐱,t)=i\mathrm{}{\displaystyle \frac{}{t}}\mathrm{\Psi }(𝐱,t)`$ (1)
where V typically confines the particle in a locale and $`V_{ext}`$ is the control taken here as $`\mu (𝐱)E(t)`$ with E(t) being the electric field and $`\mu (𝐱)`$ the dipole moment. Substituting $`\mathrm{\Psi }(𝐱,t)=A(𝐱,t)e^{iS(𝐱,t)/\mathrm{}}`$, where A and S are real functions, into Eq.(1) and separating the real and imaginary parts one easily obtains two equations. The imaginary part yields the continuity equation

$`{\displaystyle \frac{\rho }{t}}+(\rho 𝐯)=0`$ (2)
and the real part the following equation for the phase S

$`{\displaystyle \frac{S}{t}}+{\displaystyle \frac{SS}{2m}}+V_{eff}=0`$ (3)
where $`V_{eff}=V+V_{ext}+V_q`$ with $`V_q=\frac{\mathrm{}^2}{2m}\frac{^2\rho ^{1/2}}{\rho ^{1/2}}=\frac{\mathrm{}^2}{2m}[^2ln\rho ^{1/2}+(ln\rho ^{1/2})^2]`$, $`\rho =|\mathrm{\Psi }|^2`$ and v=$`\frac{S}{m}`$.

Equation (3) has the form of the classical Hamilton-Jacobi equation with an extra ’quantum potential’ term $`V_q`$. This equation can be transformed into one for the evolution of the velocity vector v by taking the gradient to give

$`{\displaystyle \frac{}{t}}𝐯=(𝐯)𝐯{\displaystyle \frac{1}{m}}(V_{eff})`$ (4)
Defining the quantum current as

$`𝐣(𝐱,t)={\displaystyle \frac{\mathrm{}}{m}}Im[\mathrm{\Psi }^{}(𝐱,t)\mathrm{\Psi }(𝐱,t)]=\rho (𝐱,t)𝐯(𝐱,t),`$
one readily obtains the equation of motion for j by substitution of $`\frac{\rho }{t}`$ and $`\frac{}{t}`$v from Eqs.(2) and (4) as

$`{\displaystyle \frac{}{t}}𝐣=𝐯(𝐣)(𝐣)𝐯{\displaystyle \frac{\rho }{m}}V_{eff}`$ (5)
Eqs.(2) and (3) or (2) and (5) describe the motion of a quantum particle within the QFD formulation of quantum mechanics. The motion of a quantum particle is governed by the current vector j and the density $`\rho `$ in Eqs.(2) and (5). Although the QFD equations resemble those of classical fluid dynamics, their quantum identity prevails due to the presence of the potential $`V_q`$ which has no classical analogue. Equivalently, the QFD equations may be viewed as those of a “classical” fluid with a highly non-linear constitutive law prescribed by $`V_{eff}`$. Various Eulerian or Lagrangian means can be exploited to solve the QFD equations , and available fluid dynamics codes may be adopted to treat these equations . The essential simplifying feature of the QFD equations is that $`\rho `$ and $`𝐣`$ or $`\rho `$ and S are often slowly varying, which is evident from quantum dynamics calculations , thereby permitting relatively coarse gridding. Despite the non-linear nature of the QFD equations, the general smoothness of $`\rho `$ and $`𝐣`$ observed lead to significant computational savings .

## III Control Expressed Within Quantum Fluid <br>Dynamics

Quantum OCT seeks the design of an external field to fulfill a particular dynamical objective. This section will provide the working equations for OCT-QFD to design an optimal electric field that drives a quantum wave packet to a desired objective at the target time t=T. The OCT-QFD formulation could be expressed in the usual way in terms of the Schrödinger equation where QFD would only act as a solution procedure. Here we will present a general approach by writing OCT directly in terms of QFD. As an example the control of a non-rotating diatomic molecule will be used as a simple illustration of the concepts. The treatment of a fully general target expectation value $`\mathrm{\Theta }_T=<\mathrm{\Psi }(T)|\mathrm{\Theta }|\mathrm{\Psi }(T)>`$ may be considered with QFD, but here we will only treat the common case where the operator $`\mathrm{\Theta }(x)`$ is only position x dependent. Then the goal is to steer $`\mathrm{\Theta }_T`$

$`\mathrm{\Theta }_T={\displaystyle _{x_l}^{x_r}}\mathrm{\Theta }(x)\rho (x,T)𝑑x`$ (6)
as close as possible to the desired value $`\mathrm{\Theta }^d`$.

The active spatial control interval is taken as $`x_lxx_r`$ over the time $`0tT`$ that the control process occurs. We desire to minimize the cost functional $`J_{cost}=J_{target}+J_{field}`$ where $`J_{target}`$ and $`J_{field}`$ are given by

$`J_{target}={\displaystyle \frac{1}{2}}\omega _x(\mathrm{\Theta }_T\mathrm{\Theta }^d)^2andJ_{field}={\displaystyle \frac{1}{2}}\omega _e{\displaystyle _0^T}E^2(t)𝑑t`$ (7)
with $`\omega _e`$ and $`\omega _x`$ being the positive weights balancing the significance of the two terms. The second term represents the penalty due to the fluence of the external field. The minimization of $`J_{cost}`$ with respect to E(t) must be subject to the satisfaction of the equations of motion for $`\rho `$ and $`𝐣`$ in Eqs.(2) and (5). We may fulfill this constraint by introducing the unconstrained cost functional as

$`\overline{J}`$ $`=`$ $`J_{cost}{\displaystyle _0^T}{\displaystyle _{x_l}^{x_r}}\lambda _1(x,t)[{\displaystyle \frac{\rho (x,t)}{t}}+{\displaystyle \frac{j(x,t)}{x}}]𝑑x𝑑t`$ (9)
$`{\displaystyle _0^T}{\displaystyle _{x_l}^{x_r}}\lambda _2(x,t)[{\displaystyle \frac{j(x,t)}{t}}+{\displaystyle \frac{}{x}}({\displaystyle \frac{j^2}{\rho }})+{\displaystyle \frac{\rho }{m}}{\displaystyle \frac{}{x}}(V+V_q+V_{ext})]𝑑x𝑑t`$
where $`\lambda _1(x,t)`$ and $`\lambda _2(x,t)`$ are Lagrange’s multiplier functions.

An optimal solution satisfies $`\delta \overline{J}=0`$, which is assured by setting each of the functional derivatives with respect to $`\lambda _1`$, $`\lambda _2`$, $`\rho `$, j and E to zero. The first two, i.e., the functional derivatives with respect to $`\lambda _1`$ and $`\lambda _2`$ regenerate the QFD equations in Eq.(2) and (5). The three others are obtained in the forms :

$`{\displaystyle \frac{\lambda _2}{t}}+{\displaystyle \frac{}{x}}(\lambda _2v_\lambda )+S_1[\rho ,j,\lambda _2]=0`$ (10)
$`{\displaystyle \frac{\lambda _1}{t}}+{\displaystyle \frac{}{x}}(\lambda _1v_\lambda )\lambda _2{\displaystyle \frac{}{x}}(V+V_q(\lambda _2)+V_{ext})+S_2[\rho ,j,\lambda _2]=0`$ (11)
and

$`{\displaystyle \frac{\delta \overline{J}}{\delta E(t)}}={\displaystyle _{x_l}^{x_r}}\lambda _2(x,t)\rho (x,t){\displaystyle \frac{}{x}}\mu (x)𝑑x+\omega _eE(t)=0`$ (12)
where

$`S_1=2{\displaystyle \frac{j}{\rho }}{\displaystyle \frac{\lambda _2}{x}}`$ (13)
$`S_2`$ $`=`$ $`{\displaystyle \frac{\lambda _2}{m}}{\displaystyle \frac{}{x}}(V_q(\rho )V_q(\lambda _2)){\displaystyle \frac{j^2}{\rho ^2}}{\displaystyle \frac{\lambda _2}{x}}`$ (16)
$`{\displaystyle \frac{\mathrm{}^2}{4m^2\rho ^{1/2}}}{\displaystyle \frac{^2}{x^2}}[{\displaystyle \frac{1}{\rho ^{1/2}}}{\displaystyle \frac{}{x}}(\lambda _2\rho )]`$
$`+{\displaystyle \frac{\mathrm{}^2}{4m^2\rho ^{3/2}}}{\displaystyle \frac{^2}{x^2}}\rho ^{1/2}{\displaystyle \frac{}{x}}(\lambda _2\rho )`$
and

$`V_q={\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{^2\lambda _2^{1/2}}{\lambda _2^{1/2}}}={\displaystyle \frac{\mathrm{}^2}{2m}}[^2ln\lambda _2^{1/2}+(ln\lambda _2^{1/2})^2]`$ (17)
The corresponding final conditions are

$`\omega _x[\mathrm{\Theta }_T\mathrm{\Theta }^d]\mathrm{\Theta }(x)\lambda _1(x,T)=0`$ (18)
and

$`\lambda _2(x,T)=0`$ (19)
Several other constraint expressions can be obtained by using equivalent forms of the continuity and dynamical equations. The form presented above is used in the subsequent numerical calculations. An alternative form in multi-dimensions symmetric between the QFD and Lagrange multiplier functions is presented in the Appendix.

The equations (9) and (10) for $`\lambda _2`$ and $`\lambda _1`$ respectively ressemble that of $`\rho `$ and j with the only difference being the extra source terms $`S_1`$ and $`S_2`$. The source terms depend on $`\rho `$ and j. $`v_\lambda `$ in the above equations is the ’velocity’ associated with the Lagrange’s multiplier and is given as $`v_\lambda =\frac{\lambda _1}{\lambda _2}`$. There are now two different quantum potential terms, one of which is a function of $`\rho (x,t)`$ and the other is a function of $`\lambda _2(x,t)`$. In this formalism the evolution of $`\lambda _1(x,t)`$ takes place by $`V_q(\lambda _2)`$ as well as the difference of the two types of quantum potential. In obtaining the above equations we have standardly assumed no variation of either $`\rho (x,0)`$ or j(x,0). Thus, we start from the initial value of $`\rho (x,0)`$ and j(x,0) to solve Eqs.(2) and (5). Eqs.(9) and (10) can be solved for $`\lambda _2(x,t)`$ and $`\lambda _1(x,t)`$ by integrating backward from time T using $`\lambda _1(x,T)`$ and $`\lambda _2(x,T)`$ given in Eqs.(15) and (16) respectively. The equations (2), (5), (9) and (10) are non-linear thereby calling for iteration to solve(cf., the algorithm in Section 4). Finally the desired control electric field is given from Eq.(11) as

$`E(t)={\displaystyle \frac{1}{\omega _e}}{\displaystyle _{x_l}^{x_r}}\lambda _2(x,t)\rho (x,t){\displaystyle \frac{}{x}}\mu (x)𝑑x`$ (20)
## 4 Application to HCl

The OCT-QFD formulation will be applied to manipulating the vibrational motion of HCl on the ground electronic state. The initial density $`\rho (x,0)=|\mathrm{\Psi }(x)|^2`$ was obtained from solving for the vibrational state from the equation

$`{\displaystyle \frac{\mathrm{}^2}{m}}{\displaystyle \frac{^2\mathrm{\Psi }(x)}{x^2}}+V(x)\mathrm{\Psi }(x)=E\mathrm{\Psi }(x)`$ (21)
using the Fourier grid Hamiltonian method  where m is the reduced mass of the HCl molecule and $`V(x)`$ is the truncated polynomial presented by Olgilvie

$`[V(x)=\{\begin{array}{cc}\hfill a_1(2\frac{xx_e}{x+x_e})^2[1+_{i=2}^9a_i(2\frac{xx_e}{x+x_e})^{i1}]b_1for& x<4\hfill \\ \hfill A[1tanh(x4)]^{3/2}for& 4x6.5\hfill \\ \hfill 0for& x6.5\hfill \end{array}]`$ (25)
where $`x_e`$=2.4086 a.u. is the equilibrium bond length of HCl. The parameters in a.u. entering the potential function are $`a_1=0.961914`$, $`a_2=1.362999`$, $`a_3=0.86675`$, $`a_4=0.49804`$, $`a_5=0.1727`$, $`a_6=0.2687`$, $`a_7=1.977`$, $`a_8=2.78`$, $`a_9=4.89`$, $`b_1=0.169695`$ and $`A=4.85\times 10^2`$. Since $`\mathrm{\Psi }(x)`$ is a stationary real function we have zero initial flux j(x,0)=0. The initial $`\rho (x,0)`$ is nearly a Gaussian packet centered around $`x_e`$. The dipole function for HCl is given by

$`\mu (x)=c_1[g(x)+c_2g^2(x)+c_3g^3(x)]`$ (26)
where $`g(x)=1tanh(\beta (xx_d))`$ and the parameters are $`c_1=0.279`$, $`c_2=0.905`$, $`c_3=1.029`$, $`\beta =0.687`$ and $`x_d=2.555`$. The following steps were carried out for implementation of the present OCT-QFD algorithm :

* Make an initial guess for the electric field E(t), which was zero in the present calculations.
* Solve the coupled equations, viz., Eq.(2) and (5) for $`\rho (x,t)`$ and j(x,t) respectively starting from $`\rho (x,0)`$ and j(x,0). The solution was achieved here by using the Flux-corrected transport(FCT) algorithm  modified for the purpose of solving the QFD equations . In doing so, we adopt the Eulerian numerical scheme.
* Evaluate the final value for $`\lambda _1(x,T)`$ given by Eq.(15) and set $`\lambda _2(x,T)`$=0 by Eq.(16).
* Solve Eqs.(9) and (10) for $`\lambda _2(x,t)`$ and $`\lambda _1(x,t)`$, respectively, by backward propagation using the same method as in step 2. Equations (9) and (10) have source terms which depend on $`\rho (x,t)`$ and j(x,t) calculated from step 2.
* Calculate the difference between the left and right sides of Eq.(16) for use in the conjugate gradient method  and calculate $`J_{cost}`$ from Eq.(7).
* Iterate steps 2 to step 6 until acceptable convergence is met.

The spatial range of the calculation was $`0x12`$ a.u., and the time interval was $`0tT`$ with T=2000 a.u. The total number of spatial mesh points is 64 which gives $`\mathrm{\Delta }x=0.1875`$ a.u. Similarly, the total number of time steps was 2048, which corresponds to $`\mathrm{\Delta }t=0.9765`$ a.u. No special effort was made to optimize the grid points, as the purpose here is to demonstrate the QFD-OCT formulation. The weight $`\omega _e`$ in Eq.(7) was taken as $`\frac{1}{2}`$, and $`\omega _x`$ = $`1000`$. The target operator was $`\mathrm{\Theta }=x`$ and $`\mathrm{\Theta }^d=3.0a.u.`$.

Figure 1 shows the control field in atomic units. The slightly non-zero values of the field at the beginning and end could be arrested by placing additional costs if desired. This pulse excites several vibrational states(not shown here) mainly by a sequence of single quantum transitions. Figure 2 shows the average distance $`<x>`$ as a function of time. The desired control value of $`<x>`$=3.0 a.u. at T is obtained through oscillatory motion of the packet. The packet is distorted in shape(not shown) while approximately retaining its original variance during the evolution. During the optimization process the total integrated probability density remained at unity up to a deviation of $`10^5`$. The iteration algorithm takes 10 steps to achieve the results shown here at 2 CPU mins. on an IRIX Silicon Graphics Machine(Release 6.1). Within numerical precision the results were the same as obtained by solving the original Schrödinger equation.

## 5. Conclusion

This paper presents a new QFD based approach for carrying out the optimal design of control fields with an illustration for the maniputation of the HCl molecule. Our previous work  shows the typical smooth and monotonic behaviour of the fluid dynamical variables, viz., S and v as opposed to the typical oscillations in the wave functions where the hamiltonian was time independent. In the present case where the system is driven with an optimal time-dependent external field we have calculated the spatial dependence of j, $`\rho `$, S and $`\mathrm{\Psi }`$ at t=T shown in Fig.3. The fluid dynamical variables(Fig.3 curves (a), (b) and (c)) used in the present method are relatively slowly varying spatial functions compared to the wave function(Fig.3, curve(d)) which apparently enhances the efficiency and the numerical saving of the present approach to controlling dynamics.

Although the illustration was for one dimension the QFD technique is directly extendable to higher dimensions, and a QFD wave packet calculation in four dimension has already been performed . The alternating direction method can effectively be used with QFD for high dimensions. Comparison with FFT propagation has been performed for two dimensional systems , showing that QFD is capable of providing a considerable increase in efficiency(i.e., by a factor of 10 or more). Regardless of the dimension, the key advantage of OCT-QFD arises from the expected smooth nature of QFD variables. A special circumstance will arise if the control “exactly” leads to a bound state with nodes that fully separates one spatial region from another. In practice placing a lower limit on the density of the order of the machine precision overcomes such difficulties. Future studies need to explore the full capabilities of the computational savings afforded by OCT-QFD.

## Acknowledgement

BD thanks Drs.Jair Botina and Tak-San Ho for useful discussions. The authors acknowledge support from the NSF and DOD.

APPENDIX
Two additional forms for the cost functional and associated initial/final conditions

The forms here are presented for reference as an alternative QFD approach. They have the advantages of simplicity and of giving equations for the Lagrange multiplier in the same form as the dynamical equations. The formalism for deriving the quations is through the Euler equations corresponding to the minimization of

$`I={\displaystyle _V}{\displaystyle _{t=0}^T}F(f,f_t,f,^2f)𝑑t𝑑V(A.1)`$
Here V denotes the volume in coordinate space. The corresponding Euler equations and conditions on time and space are

$`{\displaystyle \frac{F}{f}}{\displaystyle \frac{(\frac{F}{f_t})}{t}}({\displaystyle \frac{F}{f}})+^2({\displaystyle \frac{F}{^2f}})=0(A.2)`$
Initial condition: $`f(𝐱,0)=f_0(𝐱)`$; Final condition: $`(\frac{F}{f_t})|_{t=T}=f_T(𝐱)`$ (A.3)
Boundary conditions on dV:
$`f(𝐱,t)=f_B(𝐱,t)`$ or $`𝐧[\frac{F}{f}(\frac{F}{f})]=0`$; (A.4)
$`𝐧(f)=g_B(𝐱,t)`$ or $`\frac{F}{^2f}=0`$ (A.5)
Starting with the continuity and energy conservation equations given in Eqs.(2) and (3) in the text, we rewrite them as

$`A_t+A{\displaystyle \frac{S}{m}}+{\displaystyle \frac{A^2S}{2m}}=0(A.6)`$
$`AS_t+A{\displaystyle \frac{SS}{2m}}+VA{\displaystyle \frac{\mathrm{}^2}{2m}}^2A=0(A.7)`$
The use of the dynamical equations above in the cost functional in Eq.(8) becomes

$`J`$ $`=`$ $`{\displaystyle \frac{1}{2}}\omega _x(\mathrm{\Theta }_T\mathrm{\Theta }^d)^2+\omega _e{\displaystyle _{t=0}^T}E^2(t)𝑑t`$
$`{\displaystyle _V}{\displaystyle _{t=0}^T}[[\lambda _1(A_t+{\displaystyle \frac{AS}{m}}+A^2S)/2m](A.8)`$
$`+[\lambda _2(AS_t+A{\displaystyle \frac{SS}{2m}}+VA\mu E(t)A{\displaystyle \frac{\mathrm{}^2}{2m}}^2A)]]dtdV`$
where $`\mathrm{\Theta }_T=_V\mathrm{\Theta }(𝐱)A^2(𝐱,T)𝑑V`$. The corresponding Euler equations are obtained from the formulas in $`(A2)`$ for arbitrary variations of A, S, $`\lambda _1`$, $`\lambda _2`$ and $`E(t)`$ as

$`A_t+𝐯A=A𝐯/2`$
$`S_t+𝐯S/2=V+\mu E(t)+{\displaystyle \frac{\mathrm{}^2}{2m}}^2A/A`$
$`\lambda _{1t}+𝐯\lambda _1=[\lambda _1{\displaystyle \frac{.𝐯}{2}}+{\displaystyle \frac{\mathrm{}^2}{2m}}\lambda _2(^2A/A^2\lambda _2)/\lambda _2)](A.9)`$
$`\lambda _{2t}+𝐯\lambda _2=[\lambda _2{\displaystyle \frac{.𝐯}{2}}+{\displaystyle \frac{1}{2m}}\lambda _1(^2A/A^2\lambda _1)/\lambda _1)]`$
$`\omega _eE(t)+{\displaystyle _V}\lambda _2\mu A𝑑V=0`$
Following the formulas given in Eq.(A3) to (A5), the corresponding initial and final conditions become

$`A(𝐱,0)=A_0(𝐱)`$; $`S(𝐱,0)=S_0(𝐱)`$;

$`\lambda _1(𝐱,T)+2\omega _x(\mathrm{\Theta }_T\mathrm{\Theta }^d)A(𝐱,T)\mathrm{\Theta }(𝐱)`$; $`\lambda _2(𝐱,T)=0`$ (A.10)

The first two formulas in A.9 are equivalent to the Schrödinger equation. They can be transformed into various QFD forms in terms of $`\rho `$, v and j as in Eqs.(2) to (5) in the main text. The third and fourth equations in A.9 are the basic equations for the Lagrange multiplier functions. They are in the same flux conservation form as the QFD equations. Indeed, the third equation multiplied by $`\lambda _1`$ can be rearranged in the form of mass conservation for $`\mathrm{\Lambda }_1=\lambda _1^2`$ as

$`\mathrm{\Lambda }_{1t}+(\mathrm{\Lambda }_1𝐯)=[{\displaystyle \frac{\mathrm{}^2}{m}}\lambda _1\lambda _2(^2A/A^2\lambda _2/\lambda _2)]`$
The above derivation also can be obtained starting with the usual Schrödinger equation and its complex conjugate. Following this approach the cost functional below assures that the external field is real

$`J`$ $`=`$ $`{\displaystyle \frac{1}{2}}\omega _x(\mathrm{\Theta }_T\mathrm{\Theta }^d)^2+\omega _e{\displaystyle _{t=0}^T}E^2(t)𝑑t`$
$`{\displaystyle _V}{\displaystyle _{t=0}^T}[\lambda ^{}[i\mathrm{\Psi }_t+{\displaystyle \frac{\mathrm{}^2}{2m}}^2\mathrm{\Psi }V\mathrm{\Psi }\mu E(t)\mathrm{\Psi }]`$
$`+\lambda [i\mathrm{\Psi }_t^{}+{\displaystyle \frac{\mathrm{}^2}{2m}}^2\mathrm{\Psi }^{}V\mathrm{\Psi }^{}\mu E(t)\mathrm{\Psi }^{}]]dtdV`$
With the substitution $`\mathrm{\Psi }=Aexp(iS)`$, the cost functional reduces to the one in Eq.(A.8) with $`\lambda =\lambda _1+i\lambda _2`$.

## References

1. S. A. Rice, Science, 258, 412 (1992)
2. D. J. Tannor and S. A. Rice, J. Chem. Phys. 83, 5013 (1985)
3. A. P. Peire, M. A. Dahleh and H. Rabitz, Phys. Rev. A 37, 4950 (1988)
4. D. J. Tannor, R. Kosloff and S. A. Rice, J. Chem. Phys., 85, 5805 (1986)
5. R. Demiralp and H. Rabitz, Phys. Rev. A 47, 809 (1993)
6. J. Botina and H. Rabitz, J. Chem. Phys. 104, 4031 (1996)
7. S. Shi and H. Rabitz, Comput. Phys. Comm., 63, 71 (1991)
8. S. Shi and H. Rabitz, Chem. Phys., 139, 185 (1989)
9. W. Zhu, J. Botina and H. Rabitz, J. Chem. Phys., 108, 1953 (1998)
10. Y. Ohtsuki, H. Kono and Y. Fujimura, J. chem. Phys., 109, 9318 (1998)
11. J. Cao and K. R. Wilson, J. Chem. Phys., 107, 1441 (1997)
12. D. J. Tannor and S. A. Rice, Adv. Chem. Phys., 70, 441 (1988)
13. R. Kosloff, S. A. Rice, P. Gaspard, S. Tersigni and D. J. Tannor, Chem. Phys., 139, 201 (1989)
14. P. Brumer and M. Shapiro, Faraday Discuss. Chem. Soc., 82, 177 (1986)
15. P. Brumer and M. Shapiro, Annu. Rev. Phys. Chem., 43, 257 (1992)
16. T. Baumert and G. Gerber, Isr. J. Chem., 34, 103 (1994)
17. H. Rabitz and S. Shi in Advances in Molecular Vibration and Collisional Dynamics, Vol. 1A, 187 (1991)
18. Judson and H. Rabitz Phys. Rev. Lett. 68, 1500 (1992)
19. A. Assion, T. Baumert, M. Bergt, T. Brixner, B. Kiefer, V. Strehle and G. Gerber, Science, 282, 919 (1998)
20. C. J. Bardeen, V. V. Yakovlev, K. R. Wilson, S. D. Carpenter, P. M. Weber and W. S. Warren, Chem. Phys. Lett., 280, 151 (1997)
21. A. Assion, T. Baumert, V. Seyfried and G. Gerber in Ultrafast Phenomena edited by J. Fujimoto, W. Zinth, P. F. Barbara and W. H. Knox, Springer, Berlin (1996)
22. J. L. Herek, A. Materny and A. H. Zewail, Chem. Phys. Lett., 228, 15 (1994)
23. V. D. Kleiman, L. Zhu, J. Allen and R. J. Gordon, J. Chem. Phys., 103, 10800 (1995)
24. G. Q. Xing, X. B. Wang, X. Huang anf R. Bersohn, J. Chem. Phys., 104, 826 (1996)
25. A. Shnitman, I. Sofer, I. Golub, A. Yogev, M. Shapiro, Z. Chen and P. Brumer, Phys. Rev. Lett., 76, 2886 (1996)
26. (a) D. Bohm, Phys. Rev., 85, 166 (1952)
(b) D. Bohm, Phys. Rev., 85, 180 (1952)
27. D. Bohm, B. J. Hiley and P. N. Kaloyerou, Phys. Rep., 144, 321 (1987)
28. Bijoy K. Dey, Attila Askar and H. Rabitz, J. Chem. Phys., 109, 8770 (1998)
29. F. S. Mayor, A. Askar and H. Rabitz, J. Chem. Phys., 3, 2423 (1999)
30. R. Löhner, K. Morgan, J. Peraire, M. Vahdari, Int. J. Num. Methods Fluid, 7, 1093 (1987)
31. G. G. Balint-Kurti, C. L. Ward and C. C. Marston, Comput. Phys. Comm., 67, 285 (1991)
32. C. C. Marston and G. G. Balint-Kurti, J. Chem. Phys., 91, 3571 (1989)
33. J. F. Olgilvie, Proc. R. Soc. Lond., A 378, 287 (1981)
34. G. G. Balint-Kurti, R. N. Dixon and C. C. Marston, J. Chem. Soc. Faraday Trans., 86, 1741 (1990)
35. J. P. Boris and D. L. Book, Methods in Comp. Phys., 16, 85 (1976)
36. W. H. Press, B. P. Flannery, S. A. Teukolsky and W. T. Vetterling, Numerical Recipes, Cambridge University, New York, (1992)


format: text

text:
# Bunsetsu Identification Using Category-Exclusive Rules

## 1 Introduction

This paper is about machine learning methods for identifying bunsetsus, which correspond to English phrasal units such as noun phrases and prepositional phrases. Since Japanese syntactic analysis is usually done after bunsetsu identification \[Uchimoto et al., 1999\], identifying bunsetsu is important for analyzing Japanese sentences. The conventional studies on bunsetsu identification111Bunsetsu identification is a problem similar to chunking \[Ramshaw and Marcus, 1995, Sang and Veenstra, 1999\] in other languages.  have used hand-made rules \[Kameda, 1995, Kurohashi, 1998\], but bunsetsu identification is not an easy task. Conventional studies used many hand-made rules developed at the cost of many man-hours. Kurohashi, for example, made 146 rules for bunsetsu identification \[Kurohashi, 1998\].

In an attempt to reduce the number of man-hours, we used machine-learning methods for bunsetsu identification. Because it was not clear which machine-learning method would be the one most appropriate for bunsetsu identification, so we tried a variety of them. In this paper we report experiments comparing four machine-learning methods (decision tree, maximum entropy, example-based, and decision list methods) and our new methods using category-exclusive rules.

## 2 Bunsetsu identification problem

We conducted experiments on the following supervised learning methods for identifying bunsetsu:

* Decision tree method
* Maximum entropy method
* Example-based method (use of similarity)
* Decision list (use of probability and frequency)
* Method 1 (use of exclusive rules)
* Method 2 (use of exclusive rules with the highest similarity).

In general, bunsetsu identification is done after morphological and before syntactic analysis. Morphological analysis corresponds to part-of-speech tagging in English. Japanese syntactic structures are usually represented by the relations between bunsetsus, which correspond to phrasal units such as a noun phrase or a prepositional phrase in English. So, bunsetsu identification is important in Japanese sentence analysis.

In this paper, we identify a bunsetsu by using information from a morphological analysis. Bunsetsu identification is treated as the task of deciding whether to insert a “$`|`$mark to indicate the partition between two bunsetsus as in Figure 1. Therefore, bunsetsu identification is done by judging whether a partition mark should be inserted between two adjacent morphemes or not. (We do not use the inserted partition mark in the following analysis in this paper for the sake of simplicity.)

Our bunsetsu identification method uses the morphological information of the two preceding and two succeeding morphemes of an analyzed space between two adjacent morphemes. We use the following morphological information:

1. Major part-of-speech (POS) category,222 Part-of-speech categories follow those of JUMAN \[Kurohashi and Nagao, 1998\].
2. Minor POS category or inflection type,
3. Semantic information (the first three-digit number of a category number as used in “BGH” \[NLRI, 1964\]),
4. Word (lexical information).

For simplicity we do not use the “Semantic information” and “Word” in either of the two outside morphemes.

Figure 2 shows the information used to judge whether or not to insert a partition mark in the space between two adjacent morphemes, “wo (obj)” and “kugiru (divide),” in the sentence “bun wo kugiru.  ((I) divide sentences).”

## 3 Bunsetsu identification process for each machine-learning method

### 3.1 Decision-tree method

In this work we used the program C4.5 \[Quinlan, 1995\] for the decision-tree learning method. The four types of information, (i) major POS, (ii) minor POS, (iii) semantic information, and (iv) word, mentioned in the previous section were also used as features with the decision-tree learning method. As shown in Figure 3, the number of features is 12 ($`2+4+4+2`$) because we do not use (iii) semantic information and (iv) word information from the two outside morphemes.

In Figure 2, for example, the value of the feature ‘the major POS of the far left morpheme’ is ‘Noun.’

### 3.2 Maximum-entropy method

The maximum-entropy method is useful with sparse data conditions and has been used by many researchers \[Berger et al., 1996, Ratnaparkhi, 1996, Ratnaparkhi, 1997, Borthwick et al., 1998, Uchimoto et al., 1999\]. In our maximum-entropy experiment we used Ristad’s system \[Ristad, 1998\]. The analysis is performed by calculating the probability of inserting or not inserting a partition mark, from the output of the system. Whichever probability is higher is selected as the desired answer.

In the maximum-entropy method, we use the same four types of morphological information, (i) major POS, (ii) minor POS, (iii) semantic information, and (iv) word, as in the decision-tree method. However, it does not consider a combination of features. Unlike the decision-tree method, as a result we had to combine features manually.

First we considered a combination of the bits of each morphological information. Because there were four types of information, the total number of combinations was $`2^41`$. Since this number is large and intractable, we considered that (i) major POS, (ii) minor POS, (iii) semantic information, and (iv) word information gradually become more specific in this order, and we combined the four types of information in the following way:

| Information A: | (i) major POS |
| --- | --- |
| Information B: | (i) major POS and (ii) minor POS |
| Information C: | (i) major POS, (ii) minor POS and |
|  | (iii) semantic information |
| Information D: | (i) major POS, (ii) minor POS, |
|  | (iii) semantic information and (iv) word |

(1)
We used only Information A and B for the two outside morphemes because we did not use semantic and word information in the same way it is used in the decision-tree method.

Next, we considered the combinations of each type of information. As shown in Figure 4, the number of combinations was 64 ($`2\times 4\times 4\times 2`$).

For data sparseness, in addition to the above combinations, we considered the cases in which first, one of the two outside morphemes was not used, secondly, neither of the two outside ones were used, and thirdly, only one of the two middle ones is used. The number of features used in the maximum-entropy method is 152, which is obtained as follows:333 When we extracted features from all of the articles on January 1, 1995 of a Kyoto University corpus (the number of spaces between morphemes was 25,814) by using this method, the number of types of features was 1,534,701.

| No. of features | $`=`$ | 2 | $`\times `$ | 4 | $`\times `$ | 4 | $`\times `$ | 2 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | $`+`$ | 2 | $`\times `$ | 4 | $`\times `$ | 4 |  |  |
|  | $`+`$ |  |  | 4 | $`\times `$ | 4 | $`\times `$ | 2 |
|  | $`+`$ |  |  | 4 | $`\times `$ | 4 |  |  |
|  | $`+`$ |  |  | 4 |  |  |  |  |
|  | $`+`$ |  |  |  |  | 4 |  |  |
|  | $`=`$ | 152 |  |  |  |  |  |  |

In Figure 2, the feature that uses Information B in the far left morpheme, Information D in the left morpheme, Information C in the right morpheme, and Information A in the far right morpheme is “Noun: Normal Noun; Particle: Case-Particle: none: wo; Verb: Normal Form: 217; Symbol”. In the maximum-entropy method we used for each space 152 features such as this one.

### 3.3 Example-based method (use of similarity)

An example-based method was proposed by Nagao \[Nagao, 1984\] in an attempt to solve problems in machine translation. To resolve a problem, it uses the most similar example. In the present work, the example-based method impartially used the same four types of information (see Eq. (1)) as used in the maximum-entropy method.

To use this method, we must define the similarity of an input to an example. We use the 152 patterns from the maximum-entropy method to establish the level of similarity. We define the similarity $`S`$ between an input and an example according to which one of these 152 levels is the matching level, as follows. (The equation reflects the importance of the two middle morphemes.)

S $`=`$ s($`m_1`$) $`\times `$ s($`m_{+1}`$) $`\times `$ 10,000 $`+`$ s($`m_2`$) $`\times `$ s($`m_{+2}`$)  (2)
Here $`m_1`$, $`m_{+1}`$, $`m_2`$, and $`m_{+2}`$ refer respectively to the left, right, far left, and far right morphemes, and s(x) is the morphological similarity of a morpheme x, which is defined as follows:

| s(x) $`=`$ | 1 (when no information of x is matched) |
| --- | --- |
|  | 2 (when Information A of x is matched) |
|  | 3 (when Information B of x is matched) |
|  | 4 (when Information C of x is matched) |
|  | 5 (when Information D of x is matched) |

(3)
Figure 5 shows an example of the levels of similarity. When a pattern matches Information A of all four morphemes, such as “Noun; Particle; Verb; Symbol”, its similarity is 40,004 ($`2\times 2\times 10,000+2\times 2`$). When a pattern matches a pattern, such as “—; Particle: Case-Particle: none: wo; —; —”, its similarity is 50,001 ($`5\times 1\times 10,000+1\times 1`$).

The example-based method extracts the example with the highest level of similarity and checks whether or not that example is marked. A partition mark is inserted in the input data only when the example is marked. When multiple examples have the same highest level of similarity, the selection of the best example is ambiguous. In this case, we count the number of marked and unmarked spaces in all of the examples and choose the larger.

### 3.4 Decision-list method (use of probability and frequency)

The decision-list method was proposed by Rivest \[Rivest, 1987\], in which the rules are not expressed as a tree structure like in the decision-tree method, but are expanded by combining all the features, and are stored in a one-dimensional list. A priority order is defined in a certain way and all of the rules are arranged in this order. The decision-list method searches for rules from the top of the list and analyzes a particular problem by using only the first applicable rule.

In this study we used in the decision-list method the same 152 types of patterns that were used in the maximum-entropy method.

To determine the priority order of the rules, we referred to Yarowsky’s method \[Yarowsky, 1994\] and Nishiokayama’s method \[Nishiokayama et al., 1998\] and used the probability and frequency of each rule as measures of this priority order. When multiple rules had the same probability, the rules were arranged in order of their frequency.

Suppose, for example, that Pattern A “Noun: Normal Noun; Particle: Case-Particle: none: wo; Verb: Normal Form: 217; Symbol: Punctuation” occurs 13 times in a learning set and that ten of the occurrences include the inserted partition mark. Suppose also that Pattern B “Noun; Particle; Verb; Symbol” occurs 123 times in a learning set and that 90 of the occurrences include the mark.

This example is recognized by the following rules:

> Pattern A $``$ Partition 76.9% (10/ 13), Freq. 23
>
> Pattern B $``$ Partition 73.2% (90/123), Freq. 123

Many similar rules were made and were then listed in order of their probabilities and, for any one probability, in order of their frequencies. This list was searched from the top and the answer was obtained by using the first applicable rule.

### 3.5 Method 1 (use of category-exclusive rules)

So far, we have described the four existing machine learning methods. In the next two sections we describe our methods.

It is reasonable to consider the 152 patterns used in three of the previous methods. Now, let us suppose that the 152 patterns from the learning set yield the statistics of Figure 6.

“Partition” means that the rule determines that a partition mark should be inserted in the input data and “non-partition” means that the rule determines that a partition mark should not be inserted.

Suppose that when we solve a hypothetical problem Patterns A to G are applicable. If we use the decision-list method, only Rule A is used, which is applied first, and this determines that a partition mark should not be inserted. For Rules B, C, and D, although the frequency of each rule is lower than that of Rule A, the sum of their frequencies of the rules is higher, so we think that it is better to use Rules B, C, and D than Rule A. Method 1 follows this idea, but we do not simply sum up the frequencies. Instead, we count the number of examples used in Rules B, C, and D and judge the category having the largest number of examples that satisfy the pattern with the highest probability to be the desired answer.

For example, suppose that in the above example the number of examples satisfying Rules B, C, and D is 65. (Because some examples overlap in multiple rules, the total number of examples is actually smaller than the total number of the frequencies of the three rules.) In this case, among the examples used by the rules having 100% probability, the number of examples of partition is 65, and the number of examples of non-partition is 34. So, we determine that the desired answer is to partition.

A rule having 100% probability is called a category-exclusive rule because all the data satisfying it belong to one category, which is either partition or non-partition. Because for any given space the number of rules used can be as large as 152, category-exclusive rules are applied often444 The ratio of the spaces analyzed by using category-exclusive rules is 99.30% (16864/16983) in Experiment 1 of Section 4. This indicates that almost all of the spaces are analyzed by category-exclusive rules.. Method 1 uses all of these category-exclusive rules, so we call it the method using category-exclusive rules.

Solving problems by using rules whose probabilities are not 100% may result in the wrong solutions. Almost all of the traditional machine learning methods solve problems by using rules whose probabilities are not 100%. By using such methods, we cannot hope to improve accuracy. If we want to improve accuracy, we must use category-exclusive rules. There are some cases, however, for which, even if we take this approach, category-exclusive rules are rarely applied. In such cases, we must add new features to the analysis to create a situation in which many category-exclusive rules can be applied.

However, it is not sufficient to use category-exclusive rules. There are many meaningless rules which happen to be category-exclusive only in a learning set. We must consider how to eliminate such meaningless rules.

### 3.6 Method 2 (using category-exclusive rules with the highest similarity)

Method 2 combines the example-based method and Method 1. That is, it combines the method using similarity and the method using category-exclusive rules in order to eliminate the meaningless category-exclusive rules mentioned in the previous section.

Method 2 also uses 152 patterns for identifying bunsetsu. These patterns are used as rules in the same way as in Method 1. Desired answers are determined by using the rule having the highest probability. When multiple rules have the same probability, Method 2 uses the value of the similarity described in the section of the example-based method and analyzes the problem with the rule having the highest similarity. When multiple rules have the same probability and similarity, the method takes the examples used by the rules having the highest probability and the highest similarity, and chooses the category with the larger number of examples as the desired answer, in the same way as in Method 1.

However, when category-exclusive rules having more than one frequency exist, the above procedure is performed after eliminating all of the category-exclusive rules having one frequency. In other words, category-exclusive rules having more than one frequency are given a higher priority than category-exclusive rules having only one frequency but having a higher similarity. This is because category-exclusive rules having only one frequency are not so reliable.

## 4 Experiments and discussion

In our experiments we used a Kyoto University text corpus \[Kurohashi and Nagao, 1997\], which is a tagged corpus made up of articles from the Mainichi newspaper. All experiments reported in this paper were performed using articles dated from January 1 to 5, 1995. We obtained the correct information on morphology and bunsetsu identification from the tagged corpus.

The following experiments were conducted to determine which supervised learning method achieves the highest accuracy rate.

* Experiment 1

Learning set: January 1, 1995

Test set: January 3, 1995
* Experiment 2

Learning set: January 4, 1995

Test set: January 5, 1995

Because we used Experiment 1 in making Method 1 and Method 2, Experiment 1 is a closed data set for Method 1 and Method 2. So, we performed Experiment 2.

The results are listed in Tables 1 to 4. We used KNP2.0b4 \[Kurohashi, 1997\] and KNP2.0b6 \[Kurohashi, 1998\], which are bunsetsu identification and syntactic analysis systems using many hand-made rules in addition to the six methods described in Section 3. Because KNP is not based on a machine learning method but many hand-made rules, in the KNP results “Learning set” and “Test set” in the tables have no meanings. In the experiment of KNP, we also uses morphological information in a corpus. The “F” in the tables indicates the F-measure, which is the harmonic mean of a recall and a precision. A recall is the fraction of correctly identified partitions out of all the partitions. A precision is the fraction of correctly identified partitions out of all the spaces which were judged to have a partition mark inserted.

Tables 1 to 4 show the following results:

* In the test set the decision-tree method was a little better than the maximum-entropy method. Although the maximum-entropy method has a weak point in that it does not learn the combinations of features, we could overcome this weakness by making almost all of the combinations of features to produce a higher accuracy rate.
* The decision-list method was better than the maximum-entropy method in this experiment.
* The example-based method obtained the highest accuracy rate among the four existing methods.
* Although Method 1, which uses the category-exclusive rule, was worse than the example-based method, it was better than the decision-list method. One reason for this was that the decision-list method chooses rules randomly when multiple rules have identical probabilities and frequencies.
* Method 2, which uses the category-exclusive rule with the highest similarity, achieved the highest accuracy rate among the supervised learning methods.
* The example-based method, the decision-list method, Method 1 and Method 2 obtained accuracy rates of about 100% for the learning set. This indicates that these methods are especially strong for learning sets.
* The two methods using similarity (example-based method and Method 2) were always better than the other methods, indicating that the use of similarity is effective if we can define it appropriately.
* We carried out experiments by using KNP, a system that uses many hand-made rules. The F-measure of KNP was highest in the test set.
* We used two versions of KNP, KNP 2.0b4 and KNP 2.0b6. The latter was much better than the former, indicating that the improvements made by hand are effective. But, the maintenance of rules by hand has a limit, so the improvements made by hand are not always effective.

The above experiments indicate that Method 2 is best among the machine learning methods555In these experiments, the differences were very small. But, we think that the differences are significant to some extent because we performed Experiment 1 and Experiment 2, the data we used are a large corpus containing about a few ten thousand morphemes and tagged objectively in advance, and the difference of about 0.1% is large in the precisions of 99%..

In Table 5 we show some cases which were partitioned incorrectly with KNP but correctly with Method 2. A partition with “NEED” indicates that KNP missed inserting the partition mark, and a partition with “WRONG” indicates that KNP inserted the partition mark incorrectly. In the test set of Experiment 1, the F-measure of KNP2.0b6 was 99.66%. The F-measure increases to 99.83%, under the assumption that when KNP2.0b6 or Method 2 is correct, the answer is correct. Although the accuracy rate for KNP2.0b6 was high, there were some cases in which KNP partitioned incorrectly and Method 2 partitioned correctly. A combination of Method 2 with KNP2.0b6 may be able to improve the F-measure.

The only previous research resolving bunsetsu identification by machine learning methods, is the work by Zhang \[Zhang and Ozeki, 1998\]. The decision-tree method was used in this work. But this work used only a small number of information for bunsetsu identification666This work used only the POS information of the two morphemes of an analyzed space.  and did not achieve high accuracy rates. (The recall rate was 97.6%(=2502/(2502+62)), the precision rate was 92.4%(=2502/(2502+205)), and F-measure was 94.2%.)

## 5 Conclusion

To solve the problem of accurate bunsetsu identification, we carried out experiments comparing four existing machine-learning methods (decision-tree method, maximum-entropy method, example-based method and decision-list method). We obtained the following order of accuracy in bunsetsu identification.

Example-Based $`>`$ Decision List $`>`$

Maximum Entropy $`>`$ Decision Tree

We also described a new method which uses category-exclusive rules with the highest similarity. This method performed better than the other learning methods in our experiments.


text:
# References

Three-Pion Interferometry Results from Central Pb+Pb Collisions at 158 $`A`$ GeV/c

M.M. Aggarwal,1 A. Agnihotri,2 Z. Ahammed,3 A.L.S. Angelis,4 V. Antonenko,5 V. Arefiev,6 V. Astakhov,6 V. Avdeitchikov,6 T.C. Awes,7 P.V.K.S. Baba,8 S.K. Badyal,8 C. Barlag,9 S. Bathe,9 B. Batiounia,6 T. Bernier,10 K.B. Bhalla,2 V.S. Bhatia,1 C. Blume,9 R. Bock,11 E.-M. Bohne,9 Z. Böröcz,9 D. Bucher,9 A. Buijs,12 H. Büsching,9 L. Carlen,13 V. Chalyshev,6 S. Chattopadhyay,3 R. Cherbatchev,5 T. Chujo,14 A. Claussen,9 A.C. Das,3 M.P. Decowski,18 H. Delagrange,10 V. Djordjadze,6 P. Donni,4 I. Doubovik,5 S. Dutt,8 M.R. Dutta Majumdar,3 K. El Chenawi,13 S. Eliseev,15 K. Enosawa,14 P. Foka,4 S. Fokin,5 M.S. Ganti,3 S. Garpman,13 O. Gavrishchuk,6 F.J.M. Geurts,12 T.K. Ghosh,16 R. Glasow,9 S. K.Gupta,2 B. Guskov,6 H. Å.Gustafsson,13 H. H.Gutbrod,10 R. Higuchi,14 I. Hrivnacova,15 M. Ippolitov,5 H. Kalechofsky,4 R. Kamermans,12 K.-H. Kampert,9 K. Karadjev,5 K. Karpio,17 S. Kato,14 S. Kees,9 C. Klein-Bösing,9 S. Knoche,9 B. W. Kolb,11 I. Kosarev,6 I. Koutcheryaev,5 T. Krümpel,9 A. Kugler,15 P. Kulinich,18 M. Kurata,14 K. Kurita,14 N. Kuzmin,6 I. Langbein,11 A. Lebedev,5 Y.Y. Lee,11 H. Löhner,16 L. Luquin,10 D.P. Mahapatra,19 V. Manko,5 M. Martin,4 G. Martínez,10 A. Maximov,6 G. Mgebrichvili,5 Y. Miake,14 Md.F. Mir,8 G.C. Mishra,19 Y. Miyamoto,14 B. Mohanty,19 M.-J. Mora,10 D. Morrison,20 D. S. Mukhopadhyay,3 H. Naef,4 B. K. Nandi,19 S. K. Nayak,10 T. K. Nayak,3 S. Neumaier,11 A. Nianine,5 V. Nikitine,6 S. Nikolaev,5 P. Nilsson,13 S. Nishimura,14 P. Nomokonov,6 J. Nystrand,13 F.E. Obenshain,20 A. Oskarsson,13 I. Otterlund,13 M. Pachr,15 S. Pavliouk,6 T. Peitzmann,9 V. Petracek,15 W. Pinganaud,10 F. Plasil,7 U. von Poblotzki,9 M.L. Purschke,11 J. Rak,15 R. Raniwala,2 S. Raniwala,2 V.S. Ramamurthy,19 N.K. Rao,8 F. Retiere,10 K. Reygers,9 G. Roland,18 L. Rosselet,4 I. Roufanov,6 C. Roy,10 J.M. Rubio,4 H. Sako,14 S.S. Sambyal,8 R. Santo,9 S. Sato,14 H. Schlagheck,9 H.-R. Schmidt,11 Y. Schutz,10 G. Shabratova,6 T.H. Shah,8 I. Sibiriak,5 T. Siemiarczuk,17 D. Silvermyr,13 B.C. Sinha,3 N. Slavine,6 K. Söderström,13 N. Solomey,4 S.P. Sørensen,7,20 P. Stankus,7 G. Stefanek,17 P. Steinberg,18 E. Stenlund,13 D. Stüken,9 M. Sumbera,15 T. Svensson,13 M.D. Trivedi,3 A. Tsvetkov,5 L. Tykarski,17 J. Urbahn,11 E.C.v.d. Pijll,12 N.v. Eijndhoven,12 G.J.v. Nieuwenhuizen,18 A. Vinogradov,5 Y.P. Viyogi,3 A. Vodopianov,6 S. Vörös,4 B. Wysłouch,18 K. Yagi,14 Y. Yokota,14 G.R. Young7

(WA98 collaboration)

1 University of Panjab, Chandigarh 160014, India
2 University of Rajasthan, Jaipur 302004, Rajasthan, India
3 Variable Energy Cyclotron Centre, Calcutta 700 064, India
4 University of Geneva, CH-1211 Geneva 4,Switzerland
5 RRC Kurchatov Institute, RU-123182 Moscow, Russia
6 Joint Institute for Nuclear Research, RU-141980 Dubna, Russia
7 Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831-6372, USA
8 University of Jammu, Jammu 180001, India
9 University of Münster, D-48149 Münster, Germany
10 SUBATECH, Ecole des Mines, Nantes, France
11 Gesellschaft für Schwerionenforschung (GSI), D-64220 Darmstadt, Germany
12 Universiteit Utrecht/NIKHEF, NL-3508 TA Utrecht, The Netherlands
13 Lund University, SE-221 00 Lund, Sweden
14 University of Tsukuba, Ibaraki 305, Japan
15 Nuclear Physics Institute, CZ-250 68 Rez, Czech Rep.
16 KVI, University of Groningen, NL-9747 AA Groningen, The Netherlands
17 Institute for Nuclear Studies, 00-681 Warsaw, Poland
18 MIT, Cambridge, MA 02139, USA
19 Institute of Physics, 751-005 Bhubaneswar, India
20 University of Tennessee, Knoxville, Tennessee 37966, USA

## Abstract

Three-particle correlations have been measured for identified $`\pi ^{}`$ from central 158 A GeV Pb+Pb collisions by the WA98 experiment at CERN. A substantial contribution of the genuine three-body correlation has been found as expected for a mainly chaotic and symmetric source.

In nuclear and particle physics, the study of Bose-Einstein correlations between identical particles is widely used. It is an essential tool to obtain information on the size and time evolution of expanding systems created in heavy ion collisions . In particular, two-particle Bose-Einstein interferometry has been used for a sophisticated analysis of the dynamical evolution of the freeze out volume via selection on the transverse momentum and rapidity of the correlated particles.

Three-particle interferometry measurements can provide additional information on the space-time emission which is not accessible by two-particle interferometry. In particular, if the emission is fully chaotic, the three-particle interference study gives access to the phase of the source function’s Fourier transform, which is affected by the emission asymmetry. These asymmetries may be induced by source geometry, flow or resonance decays. If the source is not completely chaotic, as is likely to be the case, the interpretation is more difficult. Nevertheless, a comparison of the three-particle to the two-particle result allows to extract the relative strength of the true three-body correlation and can in principle remove obscuring effects, such as background from misidentified tracks or long-lived decay products, and therefore provide information more directly about the degree of coherence of the emission source. In this letter we present first results from three-particle interferometry in central 208Pb+208Pb collisions at the CERN SPS.

The fixed target experiment WA98  combined large acceptance photon detectors with a two arm charged particle tracking spectrometer. The incident 158 $`A`$ GeV Pb beam impinged on a Pb target near the entrance of a large dipole magnet. The results presented here have been obtained from an analysis of the 1995 data set. These data were taken with the most central triggers corresponding to about 10% of the minimum bias cross section of 6190 mb. The $`\pi ^{}`$ measurements were obtained with data from the negative particle tracking arm of the spectrometer. This tracking arm consisted of six multistep avalanche chambers with optical readout . A time of flight detector with a time resolution better than 120 ps allowed for particle identification. The rapidity acceptance ranged from $`y=2.1`$ to 3.1 with an average at 2.7, close to mid-rapidity which was 2.9. The momentum resolution of the spectrometer was $`\mathrm{\Delta }p/p=0.005`$ at $`p=1.5`$ GeV/c. Severe track quality cuts were applied, resulting in a final sample of 4.2$`\times 10^6`$ $`\pi ^{}`$, providing 7.2$`\times 10^6`$ pairs and 8.2$`\times 10^6`$ triplets. The $`\pi ^{}\pi ^{}`$ correlation analysis has been reported elsewhere .

The most common use of hadron interferometry concerns the study of pairs of identical particles. The two-particle correlation function can be defined as

$$C_2(𝐩_1,𝐩_2)=\frac{d^2N(𝐩_1,𝐩_2)/d𝐩_1d𝐩_2}{dN(𝐩_1)/d𝐩_1dN(𝐩_2)/d𝐩_2}$$

where $`𝐩_1`$ and $`𝐩_2`$ are the 3–momenta of the correlated particles. The product of single particle distributions in the denominator is usually obtained experimentally by a mixed event technique whereas the pair distribution in the numerator is constructed from all pair combinations of identical particles found in each event. $`C_2`$ is normalized to unity far away from the interference region. In the plane wave approximation for chaotic sources of identical particles $`C_2`$ can be written as

$$C_2(𝐩_1,𝐩_2)=1\pm |F_{12}|^2$$
(1)
with the + (-) sign for bosons (fermions). $`F_{12}`$ is the Fourier transform of the space-time source function $`S(x,k_{12})`$

$$|F_{12}|^2=\frac{|d^4xS(x,k_{12})\mathrm{exp}[iq_{12}x]|^2}{|d^4xS(x,k_{12})|^2}$$
(2)
with $`q_{12}=p_1p_2`$, the 4–momentum difference of the two particles, and $`k_{12}=(p_1+p_2)/2`$. $`F_{12}`$ is unity as $`q_{12}0`$. Thus the production of identical bosons is enhanced for pairs created close together in phase space with a small momentum difference $`Q_{12}\sqrt{q_{12}^2}`$.

Typically $`\pi \pi `$ correlation data are fit with a Gaussian form for $`|F_{12}|^2`$

$$C_2=1+\lambda \mathrm{exp}[Q_{12}^2R^2]$$
(3)
or an exponential form

$$C_2=1+\lambda \mathrm{exp}[2Q_{12}R]$$
(4)
where the parameter $`\lambda `$ is inserted to take into account the possibility that the source may not be fully chaotic and also that any wrongly reconstructed tracks, or tracks coming from decays of long-lived resonances, will dilute the correlations in the data.

The measurement of $`C_2`$ gives access to the radius $`R`$ of the source, but not to the phase $`\varphi _{12}`$ contained in $`F_{12}|F_{12}|\mathrm{exp}[i\varphi _{12}]`$ since $`C_2`$ is only a function of the square of the Fourier transform of the source distribution $`S`$. By contrast, the three-boson interference produced by a fully chaotic source is sensitive to the phase information of the Fourier transform of the source emission function. Indeed, for a chaotic source the three-body correlation function $`C_3`$, which is

$$C_3(𝐩_1,𝐩_2,𝐩_3)=\frac{d^3N(𝐩_1,𝐩_2,𝐩_3)/d𝐩_1d𝐩_2d𝐩_3}{dN(𝐩_1)/d𝐩_1dN(𝐩_2)/d𝐩_2dN(𝐩_3)/d𝐩_3}$$

can be written  as

$$C_3=1+|F_{12}|^2+|F_{23}|^2+|F_{31}|^2+2\text{Re}\{F_{12}F_{23}F_{31}\}$$
(5)
where $`F_{ij}`$ is the Fourier transform for the pair $`ij`$ contained in the triplet $`123`$. The genuine three-body correlation in $`C_3`$ is the term $`2\text{Re}\{F_{12}F_{23}F_{31}\}`$.

With $`F_{ij}|F_{ij}|\mathrm{exp}[i\varphi _{ij}]`$ and $`W\mathrm{cos}(\varphi _{12}+\varphi _{23}+\varphi _{31})`$ one may rewrite

$$2\text{Re}\{F_{12}F_{23}F_{31}\}=2|F_{12}||F_{23}||F_{31}|W.$$
(6)
Having determined $`|F_{ij}|`$ from the pair correlation function $`C_2`$, the measurement of $`C_3`$ gives direct access to $`W`$, the cosine of the sum of the three phases of the Fourier transforms, and hence provides complementary information on the shape of the source. Indeed, $`W`$ is related to the odd space-time moments of the source which generate the phases $`\varphi _{ij}`$ . For example, if the source is fully chaotic and symmetric, $`F_{ij}`$ are real, $`\varphi _{12}=\varphi _{23}=\varphi _{31}=0`$, and $`W`$ is equal to 1. $`C_3`$ is then fully determined by $`C_2`$, and the genuine three-particle correlation is maximum. If the source is not fully chaotic, more complicated expressions should be used, mixing the effects of fully chaotic and coherent sources . Nevertheless, Eq. 6 is valid in general with $`W`$ interpreted as a factor expressing the relative strength of the true three-body correlation compared to that expected for a fully symmetric chaotic source. Consequently, measuring $`W`$ values less than 1 does not allow to differentiate between asymmetries or coherence in the source.

Assuming a Gaussian form of the source function, which for two bosons leads to a correlation function described by Eq. 3, one obtains

$$C_3=1+\lambda \underset{ij=12,23,31}{}\mathrm{exp}[Q_{ij}^2R^2]+2\lambda ^{3/2}\mathrm{exp}[Q_3^2R^2/2]W$$
(7)
with $`Q_3^2Q_{12}^2+Q_{23}^2+Q_{31}^2`$. If the exponential form of the Fourier transform is assumed instead (Eq. 4), one obtains

$$C_3=1+\lambda \underset{ij=12,23,31}{}\mathrm{exp}[2Q_{ij}R]+2\lambda ^{3/2}\mathrm{exp}[(Q_{12}+Q_{23}+Q_{31})R]W$$
(8)
Inserting the values of $`\lambda `$ and $`R`$ obtained from the two-particle analysis into Eqs. 7 or 8, one can extract the three-particle strength information $`W`$.

The one-dimensional correlation function $`C_2`$, analyzed in terms of $`Q_{12}`$, is shown in Fig. 1. The two-track resolution of the spectrometer (2 cm) is dealt with by application of a proximity cut to each track pair. The data are corrected for the Coulomb and stong final state interactions in an iterative way , taking into account the source size obtained in the fit. The Gamow correction was abandoned as we found that it overcorrects the data for $`Q_{12}`$ in the range of 0.1 to 0.3 GeV/c. The effect of the experimental resolution, which is estimated by a full Monte-Carlo, can be approximated in the interference region by a Gaussian of constant $`\sigma `$ of 7 MeV/c both for $`Q_{12}`$ and $`Q_3`$. It is taken into account by replacing the function $`C_2(Q_{12})`$ used to fit the non-corrected data by

$$C_2^{rc}(Q_{12})=r(Q_{12},Q_{12}^{})C_2(Q_{12}^{})𝑑Q_{12}^{},$$

which is the convolution of $`C_2(Q_{12})`$ with the resolution function $`r(Q_{12},Q_{12}^{})`$ of the spectrometer. For display purposes, Fig. 1 is obtained by multiplying each data point by $`C_2^{rc}(Q_{12})/C_2(Q_{12})`$. Fitting the corrected data with $`C_2`$ gives the same results as fitting the non-corrected data with $`C_2^{rc}`$. The correlation function $`C_2`$ is seen to be non-Gaussian. Instead, it is better represented by an exponential form. For the $`C_3`$ analysis, an accurate description of the shape of $`C_2`$ is essential since the estimate of the $`W`$ factor extracted from $`C_3`$ depends on it. The exponential fit (Eq. 4) yielding $`R=7.29\pm 0.11`$ fm and $`\lambda =0.753\pm 0.013`$ is shown in Fig. 1.

Fig. 2 shows the three-pion correlation as a function of $`Q_3`$, after correction for resolution. A very strong $`\pi ^{}\pi ^{}\pi ^{}`$ correlation is observed, which is robust under all tracking criteria. The result shown is obtained for the same sample and the same cuts applied on the three pair combinations contained in each triplet as used for the two-pion interference analysis. For the $`C_3`$ result, the Coulomb correction applied to a particular triplet is the product of the Coulomb corrections used for the three pair combinations contained in that triplet. The resolution is taken into account using the same procedure as in the two-pion analysis. The resolution has a tiny effect compared to the Coulomb correction and it has been checked that the results are not affected by the order in which these corrections are applied to the data. The dashed line is a fit to a double exponential function

$$C_3=1+\lambda _1\mathrm{exp}[2Q_3R_1]+\lambda _2\mathrm{exp}[2Q_3R_2]$$
(9)
with fitted parameters $`R_1=5.01\pm 0.38`$ fm, $`\lambda _1=2.79\pm 0.32`$, $`R_2=1.72\pm 0.12`$ fm, $`\lambda _2=0.343\pm 0.072`$ and $`\chi ^2/d.o.f.=0.88`$. The three-pion correlation data cannot be well fitted by a Gaussian or a single exponential as a function of $`Q_3`$. Such non-Gaussian behaviour has been predicted by a final-state rescattering model .

In Fig. 2, the three-pion correlation data are compared to an estimate using Eq. 8 with $`W=1`$ (upper line) and $`W=0`$ (lower line). This estimate is made using triplets from mixed events with the $`\lambda `$ and $`R`$ parameters extracted from the two-pion interferometry analysis. Although the calculated contribution to $`C_3`$ from the genuine three-pion correlation is rather small and becomes insignificant for $`Q_3>6080`$ MeV/c, the experimental results clearly indicate a $`W`$ factor which lies between 0 and 1.

As proposed in Refs. , a method to extract the experimental value of $`W`$ as a function of $`Q_3`$ is to invert Eqs. 5 and 6 and rewrite $`|F_{ij}|`$ in terms of $`C_2`$ using Eq. 1 to obtain

$$W=\frac{\{C_3(Q_3)1\}\{C_2(Q_{12})1\}\{C_2(Q_{23})1\}\{C_2(Q_{31})1\}}{2\sqrt{\{C_2(Q_{12})1\}\{C_2(Q_{23})1\}\{C_2(Q_{31})1\}}}$$
(10)
In this method, the analysis must be performed in two steps. First, the $`\lambda `$ and $`R`$ parameters are determined both for the two-pion and three-pion correlations with a fit to the data of Eqs. 4 and 9 as previously explained. Then the data are analyzed again, and, for each triplet found, characterized by the value $`Q_3`$, $`W`$ is determined using Eqs. 4, 9, and 10 with the values $`Q_{12}`$, $`Q_{23}`$, and $`Q_{31}`$ corresponding to the three pair combinations contained in the triplet. For each bin in $`Q_3`$ containing $`N`$ triplets, the statistical error on the mean value $`W`$ is $`\sigma /\sqrt{N}`$, where $`\sigma ^2`$ is the variance of the $`W`$ distribution in this particular bin. The estimate of systematic errors is done by varying the different analysis cuts in the two and three-pion interference studies. It includes in particular the cuts used to identify the pions with the time of flight system. The systematic error on the Coulomb correction due to the error on the determination of the R parameter in the two-pion fit, as well as a possible 10% systematic error on the evaluation of the $`Q_{12}`$ and $`Q_3`$ resolution are also taken into account. The effects on $`W`$ of the statistical errors in the determination of $`C_2`$ and $`C_3`$ are treated as systematic errors by changing $`C_2`$ and $`C_3`$ respectively by $`\pm \sigma _{C_2}`$ and $`\pm \sigma _{C_3}`$. All of these variations are then added in quadrature.

Fig. 3 shows the $`W`$ values obtained for five bins of 10 MeV/c in $`Q_3`$. The error bars are the sum of statistical and systematic errors. The statistical errors (not shown separately in Fig. 3) are nearly negligible. The slight $`Q_3`$ dependence observed is not significant in view of the errors. The genuine three-body correlation is substantial with a weighted mean $`W=0.606\pm 0.005\pm 0.179`$ for $`Q_3<60`$ MeV/c. 111 The weighted systematic error is obtained by calculating the weighted average over the five $`Q_3`$ bins separately for each kind of systematic error. These errors are then added in quadrature. On the contrary, adding quadratically the systematic errors of the five $`Q_3`$ bins as done for weighted statistical errors, or as done in , would give $`\pm 0.097`$ instead of $`\pm 0.179`$ for the systematic uncertainty. This result should be compared to the lower $`\pi ^+\pi ^+\pi ^+`$ result of $`W=0.20\pm 0.02\pm 0.19`$ observed by the NA44 collaboration  in S+Pb minimum bias collisions at 200 GeV per nucleon.

In conclusion, we have studied the $`\pi ^{}\pi ^{}\pi ^{}`$ interference of pions produced in central Pb+Pb collisions at 158 GeV per nucleon. Although its contribution is small, the genuine three-pion correlation, the portion of the correlation not trivially due to the two-pion interference, has been extracted and found to be substantial. The genuine three-pion correlation is greater than reported for S+Pb minimum bias collisions , but not as large as expected for a fully chaotic and symmetric source.

We would like to thank the CERN-SPS accelerator crew for providing an excellent Pb beam. This work was supported jointly by the German BMBF and DFG, the U.S. DOE, the Swedish NFR, the Dutch Stichting FOM, the Swiss National Fund, the Humboldt Foundation, the Stiftung für deutsch-polnische Zusammenarbeit, the Department of Atomic Energy, the Department of Science and Technology and the University Grants Commission of the Government of India, the Indo-FRG Exchange Programme, the PPE division of CERN, the INTAS under contract INTAS-97-0158, the Polish KBN under the grant 2P03B16815, the Grant-in-Aid for Scientific Research (Specially Promoted Research & International Scientific Research) of the Ministry of Education, Science, Sports and Culture, JSPS Research Fellowships for Young Scientists, the University of Tsukuba Special Research Projects, and ORISE. ORNL is managed by UT-Battelle, LLC, for the U.S. Department of Energy under contract DE-AC05-00OR22725.