format: text

text:
# SPINning Parallel Systems Software This work was supported by the Mathematical, Information, and Computational Sciences Division subprogram of the Office of Advanced Scientific Computing Research, U.S. Department of Energy, under Contract W-31-109-Eng-38.

## 1 Introduction

Reasoning about parallel programs is surprisingly difficult. Even small parallel programs are difficult to write correctly, and an incorrect parallel program is equally difficult to debug. In our experience writing the Multi Purpose Daemon (MPD), a parallel system program described below, this characterization has been borne out: despite MPD’s small size and apparent simplicity, errors have impeded progress toward code in which we have complete confidence. (Of course, another possible explanation is that as parallel programmers we are inept; we reject this hypothesis because it is impossible to verify rigorously.)

Such a situation motivates us to explore program verification techniques. Since our programs are small and our algorithms simple (when viewed from the perspective of a single process), we hope that program verification software will be able to handle our problem; since the difficulty of reasoning about parallelism has shown us that we really do need help with this problem, investing in verification is worth the effort.

MPD \[bgl00:mpd:pvmmpi00, butler-lusk-gropp:mpd-parcomp\] is a process manager for parallel programs and is itself a parallel program. Its function is to start the processes of a parallel job in a scalable way, manage input and output, deal with faults, and cause jobs to terminate cleanly. While the job is running, it may need to provide services to the application, such as implementing a barrier or assisting an application process in setting up communication with another process in the job. MPD is the sort of process manager needed to run applications that use the standard MPI \[mpi-forum:journal, mpi-forum:mpi2-journal\] library for parallelism, although it is not MPI-specific. MPD is distributed as part of the portable and publicly available MPICH \[gropp-lusk:mpich-www, gropp-lusk-doss-skjellum:mpich\] implementation of MPI.

Our first attempt \[mpd-acl2\] to use formal verification techniques to ensure correctness of MPD algorithms was based on the ACL2 \[acl2:book, acl2\] theorem prover. While ACL2 provided a useful simulation environment, formulating desired properties of and reasoning about models of MPD algorithms proved difficult. Our second approach, described here, employs the model checker Spin \[spin:book, spin:article\]. Our particular application is unusual in that the number of entities and the topology of the communication network can change over time. We also need to model a larger number of entities than many Spin applications do. In Section 4 we discuss the challenges that these properties provide to the Spin system.

In Section 2 we describe MPD in more detail and outline our method for modeling a distributed, dynamic set of Unix processes in Promela. In Section 3 we describe our experiences with this approach, which we believe shows potential benefits for the further development of MPD. In Section 4 we present the concrete results of specific verification experiments, and we conclude in Section 5 with a summary of the current project status and our future plans.

## 2 Approach

To present what we want to verify and how we have gone about it, we describe here a few salient features of the MPD system. More details can be found in \[bgl00:mpd:pvmmpi00\] and \[butler-lusk-gropp:mpd-parcomp\].

The MPD system consists of several types of processes. The daemons are persistent (may run for weeks or months at a time, starting many jobs), and there typically exists one daemon instance per host in a TCP-connected network. Manager processes are started by the daemons to control the application processes (clients) of a single parallel job and provide most of the MPD features. The daemons are connected in a ring. A console process is started by a user or another process to connect to the daemon ring and give it a command, such as mpirun to start a user parallel job. Separate managers for each user process, started by the daemons, support individual process environments for the user processes. The managers also connect themselves into a ring.

Exactly how the daemons are started or connected is not important, since the system provides a number of choices, and the process need not be particularly fast. A console command is started by the user, either interactively or under the control of a batch scheduler. The daemons spawn the managers, which use information given them by the daemons to connect themselves into a ring, then spawn the clients. The startup messages traverse the ring quickly, so most invocation of new processes and connecting takes place in parallel, leading to fast startup even for jobs involving hundreds of processes.

The situation is as shown in Figure 1, where the clients may be application MPI processes. The vertical solid lines represent connections based on pipes; the remaining solid lines all represent connections based on Unix sockets. The dashed lines represent the trees of connections for forwarding input and output, and the dotted lines represent potential connections among the client processes. The dot-dashed line is the original connection from console to local daemon on a Unix socket, which is replaced during startup by the network connection to the first manager.

An important feature of MPD for our purpose here is that the structure of the code for each of the three process types (daemon, manager, console) is essentially the same: after initialization, the process enters an essentially infinite loop, implemented by a Unix socket function select, which indicates which sockets have messages available. That is, most of the time, the process is idle. When a message arrives on one of its sockets from one of the processes it is connected to, it wakes up, parses the message, calls the appropriate handler routine to process the message, and re-enters the idle state by calling select again. The handler routine itself does a small amount of processing, typically resulting in the creation of new sockets or in the sending of messages on existing sockets. The logic of the distributed algorithms executed by the system as a whole, therefore, is contained primarily in the handlers. No individual message, received by one process, results in much activity. The individual handlers are typically implemented in a few lines, or at most a few tens of lines, of C code.

This structure allows us to treat the system as comprising three layers (see Figure 2). The top tier corresponds to the upper-level logic of the process (initialization, select loop, parsing and dispatching of incoming messages to handlers). This is sequential logic that we have confidence in; hence, although we must model this part of each process in some way in Promela, the Promela model does not have to be faithful to the algorithms in this layer. The bottom layer corresponds to well-understood Unix operations on sockets. Again the code is sequential and not of particular interest. It is in the middle layer (the handlers) that the interesting parallel algorithms are expressed and the bugs appear.

Verification of the algorithms executed by different components of the MPD system is not the only goal of the project. Equally important and desirable is the ability to translate, possibly automatically, a verified Promela model to executable C code while preserving the verified properties of the model. Only the middle layer of the model has to be translated. Thus, in our Promela model of each type of process, we take considerable liberties with the top and bottom layers, but we wish the Promela code for the middle layer (the handlers) to be directly translatable into C or the scripting language Python. Our current implementation of MPD is in C.

### 2.1 Modeling Components of the MPD System

Components of the MPD system map naturally to Promela entities: a proctype is defined for each different MPD process type; individual daemon, manager, console, and client processes correspond to active instances of the corresponding proctypes; sockets map to channels; and messages that are read and written over the sockets correspond to messages traveling on the channels. Our models of individual process types preserve the three-layer structure (Figure 2) of executable MPD code for clarity, readability, and modularity and also to facilitate the translation of verified handler algorithms to executable code.

The top tier of the model typically contains the start-up logic of a particular MPD algorithm followed by a mechanism for sampling and processing the attached channels. The middle tier corresponds to a collection of handlers that define the set of message types permitted for each type of communication link as well as the behavior of the process in response to each message type. The lowest tier consists of a Promela-based implementation of the Unix socket primitives. Note that an MPD programmer relies on the predefined functionality of the socket processing functions but does not implement the primitives. Creating a Promela library of Unix socket primitives allows us to (1) hide the details of the socket model from both the verification and translation to the executable code, (2) interchange, if need be, different models of sockets without changing the remainder of the model, and (3) reuse the socket model in verification of independent MPD algorithms.

### 2.2 Modeling Unix Sockets

A Unix socket is an endpoint of a bidirectional communication path between two processes. In MPD, sockets are manipulated when a connection between two processes is established or destroyed. Correct operation of some MPD algorithms depends on correct allocation and manipulation of sockets, while other MPD algorithms assume a static system of processes and communication links between them. Our objective is to model sockets efficiently to support verification of both types of MPD algorithms.

An MPD process runs on top of an operating system, which provides (among other services) the implementation of sockets and their handling. We must model the sockets, but we do so without creating a model of the operating system. In fact in our model, the operating systems of the hosts on which MPD processes execute are combined into a single implicit global operating system. The only explicit manifestation of the operating system is the socket descriptor structure described below. In the model, the functionality of the operating system is hidden inside the socket library, and its tasks are handled directly by the MPD processes.

A Unix socket is referenced by a file descriptor (fd) and represents a buffer for reading and writing messages. Our Promela model of a socket consists of a channel and a three-part record that describes how the particular socket should be used. A model of an MPD algorithm contains an array of channels and an array of socket descriptor structures. The first field of the socket descriptor structure references the fd at the other endpoint of the connection. The second field identifies a process that has an exclusive privilege to read from, write to, and deallocate the socket. The third field is a flag that indicates whether the socket has been allocated to a particular process and, if so, how it can be used. Outside the context of MPD, the usage flag need indicate only whether the socket is free or allocated. However, since the model is created specifically for use with MPD, the flag field of the socket descriptor structure is also used to denote how the socket is used by an MPD process. For example, when a ring of MPD processes is established, sockets corresponding to the right-hand side of the connection are processed differently from the sockets on the left. Figure 3 shows two connected MPD processes. Figure 4 shows the corresponding state of an array of socket descriptors as well as one unallocated socket.

The following Unix socket primitives have been modeled in accordance with their defined functionality\[stevens-unp1\]: connect, accept, close, read, and write. Other socket primitives, such as socket, bind, and listen, are not explicitly modeled but are implied in the model. Select is used only in the top tier of the MPD process. We include a model of select in the socket library, but its implementation is greatly influenced by the specific way in which the function is used by MPD algorithms. The modeled socket functions, defined as inline Promela functions of the same name, serve as an interface between the bottom and middle tiers of the model. Below is an excerpt of the socket library together with the definitions of the channels and socket descriptor structure. Definitions of self-explanatory functions and macros are omitted.

```
chan connection[CONN_MAX] = [QSZ] of {msg_type};

typedef conn_info_type {
  unsigned other_fd  : FD_BITS;
  unsigned owner_pid : PROC_BITS;
  unsigned use_flag  : FLAG_BITS;
  };

conn_info_type conn_info[CONN_MAX];

inline read(file_desc, message)
{
  connection[file_desc]?message;
}

inline write(file_desc, message)
{
  connection[conn_info[file_desc].other_fd]!message;
  fd_select_check(conn_info[file_desc].other_fd)
}

inline close(file_desc)
{
   IF  /* other side has not been closed yet */
   :: (conn_info[file_desc].other_fd != INVALID_FD) ->
      set_other_side(conn_info[file_desc].other_fd,INVALID_FD);
      fd_select_check(conn_info[file_desc].other_fd)
   FI;
   deallocate_connection(file_desc)
}

inline connect(file_desc, lp)
{
   allocate_connection(j);          /* server’s connection */
   set_owner(j, lp);
   set_handler(j, AWAIT_ACCEPT);
   allocate_connection(file_desc);  /* client’s connection */
   set_owner(file_desc, _pid);
   set_other_side(j, file_desc);    /* relate connections */
   set_other_side(file_desc, j);    /* to each other      */
   lp_select_check(lp)
}

inline accept(file_desc)
{
   file_desc = 0;
   do
   :: (file_desc >= CONN_MAX) ->
      assert(0) /* block if no connect was made */
   :: (readable_lp(file_desc,_pid)) ->
      set_handler(file_desc, NEW);
      break
   :: else ->
      file_desc = file_desc + 1
   od
}

```

Notice that read(fd) and write(fd) operate on different sockets. The reading operation amounts to receiving a message from the connection\[fd\] channel, while the writing operation places the message on the other buffer of the connection, which is pointed to by the other\_fd field of conn\_info\[fd\]. This approach ensures that the connection is truly bidirectional and that the two processes at the endpoints of the connection can read and write independently. Connect and accept are companion operations. To establish a connection, a client process connects to a listening port on the server process. To complete the connection, the server must accept the connection request from the client. In the model, the client allocates both sockets of the connection and sets the use\_flag of the server’s socket to AWAIT\_ACCEPT. The accept operation then locates the socket with the set flag. In Unix, an accept will block if executed before a connect is issued. The Promela model also implements this behavior. Finally, the close operation deallocates the socket. When a close is issued in Unix, the process on the other side of the connection essentially sees an EOF on the corresponding socket. To simulate the behavior in Promela, the model of the close operation sets the other\_fd of the remaining part of the connection to an invalid value. Write, connect, and close operations contain references to {fd,lp}\_select\_check functions, which are helper functions to the select operation. They set a bit for the owner process of the connection, indicating that the process can leave the idle state and do some processing in response to received messages or other set flags. The select function blocks until the bit is set.

### 2.3 The Road Not Taken

One may model the MPD system in many ways. We are fully prepared for the possibility that as this project progresses and we model other kinds of MPD algorithms, our model will change. The three-layered approach allows us to make dramatic changes in different levels of the model without affecting the remainder of the model. We experimented with several different models of the MPD system and came to the following conclusions.

While it is tempting to model the operating system explicitly and to hide the manipulation of the sockets inside such a model, doing so goes against the methodology of Spin/Promela model construction. An explicit model of the operating system produces a process that simply forwards messages between daemons, managers, clients, and consoles, resulting in a rapid explosion of the state space.

One can model the socket operations by sending explicit messages. For example, a connect can result in a message that is consumed by the corresponding accept. A close operation can send a special eof message on the other end of the connection, if it is still open. The select operation may be viewed as a message whose parameters include references to sockets, which should be sampled for messages. The select would return when a response message is received. The parameters of the response message could contain the references to the sockets that in fact have messages that are ready for consumption. After experimenting with such an approach, we decided that the model should contain only those messages that represent explicit communication between MPD processes, while communication that occurs between the host operating systems should be represented by other means, such as setting and resetting of flags.

## 3 Early Experiences with Modeling and Verification

This is our first Spin project. Initial efforts concentrated on investigating applicability of the Spin-based approach to our problem. We experimented by constructing different models of Unix sockets and different models of MPD algorithms and attempting verification of these models. These efforts resulted in an early success in demonstrating that a proposed MPD algorithm was incorrect.111Admittedly, our focus at the time was on learning Spin and Promela and not on development of MPD algorithms. Had we been concentrating on the latter topic, the rather obvious error we discovered with Spin’s aid probably would have been found during early stages of design. Nonetheless, we view the experience as evidence of the usefulness of the Spin-based approach to verification of MPD algorithms.

### 3.1 A Buggy Algorithm for Creating a Daemon Ring

Establishment of a ring of daemons, the first step of the MPD system, and maintenance of the ring are central to the operation of MPD. Informally, daemon ring creation proceeds as follows. The initial daemon establishes a listening port to which subsequent connections are made. The daemon connects to its own listening port, creating a ring of one daemon. The listening port of the first daemon and the name of the host processor are queried from the console. The desired number of daemons is then initiated and directed to enter the ring by connecting to the first daemon. Figure 5

shows the result of inserting a new daemon into an existing ring. Upon completion of the insertion the old connection between daemons on the right and left of the new daemon is disconnected (shown in the figure by the dashed line). Note that in the special case of insertion into a ring of one daemon, the daemon plays both the left and the right roles.

Initially we modeled an algorithm that allowed new daemons to enter the ring sequentially. Figure 6

presents the algorithm as a message sequence chart. The new daemon is supplied with the information it needs to initiate a connection to the daemon D-left, so identified because it will end up as the left-hand side neighbor of the new daemon. The connection is accepted by D-left. The new daemon then queries the D-left daemon about the coordinates of its right-hand side neighbor. As no information about the state of the ring is kept by the daemons in this algorithm, D-left queries D-right for its listening port number and the name of the host. D-right returns the requested information, which is then forwarded to the new daemon. The new daemon declares itself the new right-hand side neighbor of D-left, which replaces its existing right-hand side connection with the connection to D-new. The old connection is closed. The new daemon continues to enter the ring by connecting to D-right using the coordinates it received in the rhs\_info\_return message from D-left. The connection is accepted, and the new daemon declares itself the new left-hand side of D-right, which closes the remaining endpoint of the old connection.

MPD design requirements, however, mandate that several ring insertion commands can be issued simultaneously. Moreover, an MPD process can be blocked only by a select. (The second requirement is fulfilled if there are no blocking read statements within the body of any handler.) Our sequential algorithm exhibits neither of the mandated properties. The algorithm assumes that the sequence of messages rhs\_info\_request, …, new\_rhs is not interrupted, that is, D-new sends an rhs\_info\_request message and blocks until it receives the rhs\_info\_return message. D-left is similarly blocked upon sending the rhs\_info\_return message until the receipt of the new\_rhs message. If the algorithm is changed so that daemons are blocked only by a select and several new daemons simultaneously enter the ring, an erroneous execution scenario is possible. The error occurs when two daemons attempt to simultaneously enter a ring of more than one daemon. The initial connection is made by both new daemons to the same daemon in the ring. Then both daemons issue an rhs\_info\_request. If the processing of the resulting messages is interleaved, both new processes are told to connect their right-hand side to the same daemon. A ring will not be established in this case. In the correct algorithm, one of the new processes should connect its right-hand side to the other new daemon, which should connect its right-hand side to the existing daemon. The algorithm is fixed by storing on each daemon the coordinates of its right-hand side neighbor. Doing so makes unnecessary the series of messages to request and return this information to the new daemon.

This error was discovered during simulations of the Promela model within Spin. The experience was important not only because Spin helped us find an error but also because Spin helped us to do so very quickly. This experience leads us to believe that the tool will become invaluable in future development efforts both for verification and for rapid prototyping and testing of new algorithms.

### 3.2 Early Problems with Modeling and Verification

Our early models were naive because of our inexperience and our attempts to match the model too closely to C code. Including an explicit model of the operating system for socket handling led to extreme state space explosion while requiring a long state vector size. After examining the literature on pragmatic use of Spin \[spin:abstraction, spin:low-fat\], we developed much leaner models. In optimizing our original explicit models we tried to find a level of abstraction that preserved the correlation of the Promela model to the eventual C code while keeping verification feasible.

## 4 Verification of MPD Algorithms

An MPD daemon ring is a dynamic structure: new processes may enter the ring at any time, and existing processes or their host processors may crash or may shut down in an orderly fashion. We modeled and verified algorithms for ring creation and ring recovery after a single process/processor crash. The majority of MPD algorithms reside in the managers. We modeled and verified a barrier algorithm, an example of a manager-level functionality, which ensures that all clients reach a certain point in the execution of a parallel job before any client is allowed to proceed further.

All verification runs were conducted on a 933 MHz Pentium III processor with 970 MB of usable RAM. We used default XSpin settings for all verification attempts, except when we increased the memory limit to allow the search to complete. In cases where verification did not complete with default parameters within physical memory limits, verification with compression (-DCOLLAPSE compile-time directive) was performed. Such experiments are identified by an asterisk.

### 4.1 Ring Establishment Algorithm

Figure 7 shows a message-sequence chart

representation of an algorithm that allows parallel insertion of daemons into the ring. A minor difference between this algorithm and the ring establishment algorithm described in Section 3 (Figure 6) is that the new\_rhs message is sent before the coordinates of D-right are supplied to D-new via the message reconnect\_rhs. The major difference is that in this algorithm each daemon records and maintains the identities (listening port number and the host name) of the two daemons located to its right. Strictly speaking, only the identity of the immediate right-hand side neighbor is needed to establish the ring. To recover from an unexpected crash of a daemon, however, the identity of the second right-hand side neighbor (rhs2) must be known as well. The algorithm involves sending a message of type rhs2info, to communicate this information to the appropriate daemon. This operation is not shown in the figure. In a small ring (fewer than three daemons) or when a ring is being created, special case logic ensures that appropriate information is forwarded to the correct daemon. In the general case, upon receipt of the new\_rhs message D-left sends the rhs2info message, in a counterclockwise direction along the ring, to its left-hand side neighbor, indicating that it has a new rhs2, namely, D-new. In addition, upon receipt of the new\_lhs message, D-right sends the rhs2info message along the ring, also in a counterclockwise direction, to D-new. The rhs2 neighbor of D-new is the right-hand side neighbor of D-right.

The algorithm was modeled by using the three-tiered approach described above. To check whether the algorithm is correct, we verify that the resulting system topology, which is implicit in the socket descriptor structures array, is in fact a ring of the correct size. In addition, we check that the state information (identities of the two neighbors to the right) agrees with the information in the socket descriptor structures array. These two checks are performed when the ring establishment algorithm is completed, which in reality corresponds to all daemons entering an idle state. In the Promela model this corresponds to a timeout. A properly connected ring manifests itself in successfully supporting operation of daemon-level algorithms. In fact, MPD designers and users test whether the ring was established successfully by executing one of such algorithms and examining its results. The algorithm, invoked by mpitrace, reports the identities of all daemons and their position in the ring. To convince our “customers” (i.e., MPD designers) that the Promela model of the ring establishment algorithm is correct, we also verify that upon its completion, the trace algorithm terminates after having visited every daemon in the ring. The trace algorithm is implemented both in MPD and in our model by sending two kinds of messages along the ring, which leads to an increase in the number of states that have to be examined during a verification attempt. For this reason, on a model of a given size, verification of the successful trace completion is more expensive, in terms of time and memory, than verification of ring topology and state information.

Statistics for the verification of the ring algorithm are presented in Table 1. The first portion of the table reports on the verification of state of the ring properties, that is, the information in the socket description structures array and the recorded identities of the two right-hand side neighbor of each daemon. The second portion of the table presents statistics on verification of a successful trace after ring establishment. The algorithm was verified for models comprising up to four daemons, as indicated in the second column of the table.

We were unable to exhaustively verify the algorithm on models with five or more daemons. Because of dramatic state space explosion and large state vectors, verification attempts for these models ran out of 970 MB of memory, and using the compression and graph encoding techniques (-DCOLLAPSE and -DMA=n compile-time directives) still did not enable the search to complete. When we applied predicate abstraction techniques to model the ring establishment algorithm, verification succeeded for models with only up to eight daemons. However, the desired correlation of the Promela model to the C/Python code was lost, as was the ability to perform meaningful simulations. As shown below in verification statistics of other MPD algorithms that we have modeled and verified, the rapid increase in the number of states for relatively small models occurs only for the ring establishment algorithm and is due in large part to its properties.

Since the daemons enter the ring in parallel, there are many possible interleavings of their execution leading to many possible ring configurations. In general, given $`n`$ daemons, there are $`n!`$ resulting ring configurations. Moreover, a daemon enters the ring in two steps: first it connects its left-hand side, then its right-hand side, which in turn increases the number of possible interleavings. An increase in model size also leads to an increase in the number of sockets/fds that are manipulated. Since every connection consists of two sockets, there is a minimum of $`2n`$ sockets in a ring of $`n`$ daemons. However, as the ring is being created, additional sockets are required because there are execution sequences in which allocation of a new socket occurs before deallocation of a socket belonging to a connection whose other end has been closed. Finally, an increase in model size leads to an increase of the buffer size of the communication channels. Recall that when a new daemon enters the ring, a rhs2info message is sent along the ring in a counterclockwise direction to a daemon located two positions to the left of the new daemon. Therefore, in a general case, for $`n`$ new daemons inserted into the ring, $`n`$ rhs2info messages are sent. Some execution sequences result in all $`n`$ messages accumulating on a single communication channels. Therefore, the buffer size is equal to the number of daemons.

The fact that verification succeeds only on small models needs to be put in perspective. On one hand, a running MPD typically consists of several tens or hundred of processes. Therefore, given the current technology, we cannot completely verify MPD algorithms for every possible system size and topology. On the other hand, prior experience with debugging of the MPD code suggests that even the most difficult errors manifest themselves in systems of just a few (four to ten) processes. Therefore, the models of MPD algorithms as the current level of abstraction allow us to perform verification of some algorithms on models of satisfactory size. For other algorithms, such as the ring establishment algorithm, a slightly more abstract model or a more efficient socket library will enable verification to complete on models of sufficient size.

### 4.2 Recovery from a Single Nondeterministic Failure

The ring recovery algorithm works as follows. When a daemon in a properly established ring fails, the operating system on its host processor will close all sockets that belonged to it, so the neighbors of the failed daemon will see that the sockets on their opposite end of their connection to the daemon were closed. Of course, sockets may be closed for legitimate reasons as well, but the suite of MPD algorithms is designed in such a way that a closed socket on the right-hand side connection in the ring, without any advance notification that a controlled shutdown or disconnection will take place, signifies an unintended failure and places the burden of recovery on the remaining members of the ring. In the algorithm, the left-hand side neighbor reinstates the ring by establishing a connection to the right-hand side neighbor of the failed daemon. The identities of the two right-hand side neighbors must be updated for all affected daemons.

In our model of the algorithm the initial ring is hard-coded. One daemon is directed to fail, although we do not specify which one should do so. The recovery procedure is then initiated. The model algorithm was verified against the three correctness properties of a properly established ring, as discussed above. Table 2 shows statistics of the verification attempts.

### 4.3 Manager-Level Barrier Algorithm

Parallel programs frequently rely on a barrier mechanism to ensure that all processes of the job reach a certain point (complete initialization, for example) before any are allowed to proceed further. Parallel jobs, that is, programs running on the clients, rely on the manager processes to implement the barrier service. The algorithm proceeds as follows. A manager process is designated as the leader of the algorithm and is given a rank of 0. When the leader reads a request from its client to provide the barrier service, it sends a message barrier\_in to its right-hand side neighbor in the ring. When a non-leader manager receives the barrier\_in message, its behavior is determined by whether its client has already requested the barrier service. If the client has done so, the manager forwards the message to the right-hand side manager. Otherwise, it holds the barrier\_in message until the request from the client arrives. While the barrier\_in message is held, a bit variable holding\_barrier\_in is set. Once the barrier\_in message traverses the entire manager ring and arrives back in the leader, meaning that each client has reached the barrier and notified its manager, the leader sends a barrier\_out message around the ring. When a manager receives the barrier\_out message, it notifies its client to proceed past the barrier. The leader can be either the first or the last manager to allow its client to proceed.

We modeled the algorithm on top of the socket library. A ring of managers is hard-coded in the beginning, and a manager with \_pid of 0 is designated as the leader. There is no need to model the clients explicitly; thus the communication with them is represented by two bits per manager. One bit designates that a request for a barrier service has been received by the manager. The other bit, client\_barrier\_out, designates that the manager notified the client to proceed past the barrier. Global arrays of bits were defined, using the bit-array implementation by Ruys \[spin:low-fat\], to store these values. A constant ALL\_BITS corresponds to a value of the bit array where an element was set for every manager in the ring.

We verified two correctness conditions about the algorithm. First, at the end of the algorithm, all clients must have been told to proceed past the barrier:

```
timeout -> assert(client_barrier_out == ALL_BITS)

```

The second condition is an invariant: no client is allowed to proceed until all clients have reached the barrier and all managers have released the barrier\_in message:

```
assert((client_barrier_out == 0) ||
       ((client_barrier_in == ALL_BITS)
        && (holding_barrier_in == 0)))

```

Table 3 shows statistics for verification of the barrier algorithm. We were able to exhaustively verify models with up to fourteen managers.

## 5 Summary and Future Plans

We described here our first experiences in applying the Spin-based approach to verification of a parallel process management system called MPD. We settled on a three-tier architecture for the models in order to maintain the correlation to the eventual code of the MPD system and to enforce modularity of the model. The bottom layer of the architecture consists of a Promela model of operations on Unix sockets. We encountered some early difficulties in the verification attempts. Specifically, for the ring establishment algorithm, exhaustive verification can be completed only on models with up to four daemons. However, we were able to exhaustively verify larger models of other algorithms.

Based on our experiences, we believe that design and development of algorithms for MPD and similar systems can benefit greatly from application of the Spin-based software verification methods. Spin’s simulation capability allows for rapid prototyping of new algorithms. Since even the most difficult errors can be discovered on models comprising only a few processes, the verification engine of Spin enables us to verify the algorithms on models that are sufficiently large for our purposes.

A long-term goal of this project is to model and verify MPD algorithms and then translate them into C or another programming language, while preserving the verified properties. Ideally, translation should be automated. To allow this to happen, the Promela model must not be overly abstract. Figure 8 shows a Promela

model and a C implementation of a portion of the barrier algorithm, in which a barrier\_in message is received and processed by a manager. Automated translation certainly appears feasible for this level of abstraction of the Promela model. Notice the one-to-one correspondence between the control structures of the two segments. There is further correspondence between Promela and C for checking and setting of boolean variables (lines p.3,c.2 and p.15,c.13). The code for message assembly (lines p.6,c.4 and p.10,c.12) matches as follows. A Promela macro make\_barrier\_out\_msg corresponds in C to writing a string, containing a barrier\_out command, to a buffer. In the Promela model, before a message can be written to a buffer, the corresponding fd must be found using functions find\_right or find\_left. In the C code, the buffer is written to a file descriptor correspondingly referenced by rhs\_idx or lhs\_idx. Therefore, the two Promela lines p.6-7 match a single C line c.5.

Of course, not all models will lend themselves well to verification at this level of abstraction, as verification of the ring establishment algorithm demonstrated. But, models of the MPD algorithms should fall into just a few different classes with respect to the level of abstraction, and a separate mapping can be defined for each such level to enable the Promela-to-C translation.

We will continue to model and verify individual MPD algorithms. They include the daemon-level functionality for controlled shutdown of a portion of the ring and subsequent ring reestablishment, as well as various manager-level algorithms, such as the handling of the parallel input and output to the console. Correct interaction of these algorithms is also very important. Many things take place in parallel in a running MPD system. Daemons enter and leave the ring, as do managers, different client processes request different services from the managers, and several instances of the same algorithm may be executing simultaneously. We hope to be able to reason formally about MPD models that consist of several related and interdependent algorithms.

The Promela implementation of the the Unix socket library as well as models of the MPD algorithms described in this paper are available at http://www.mcs.anl.gov/ matlin/spin-mpd.


text:
# 1 Introduction

## 1 Introduction

The technique of partial reconstruction of $`D^{}`$ mesons (charge conjugate states are always implied), in which only the slow pion from $`D^{}\overline{D}{}_{}{}^{0}\pi _{}^{}`$ is reconstructed, has been widely used in the past  to select large samples of reconstructed $`B`$ mesons. This technique provides a way to measure the combination of the CKM unitarity triangle angles $`(2\beta +\gamma )`$ with $`B^0D^{}\pi ^+`$ decays . The present measurement has been performed as a first step towards the goal of an analysis measuring the angle $`\gamma `$. In this respect, the reconstruction of the signal events, the rejection of the background, the characterization of the various background components and finally the study of the $`\mathrm{\Delta }t`$ resolution obtained in these events are important tools for a CP analysis. All these tools are presented here in the context of a $`B^0`$ lifetime measurement.

## 2 The BABAR detector and dataset

The data used in this analysis were collected with the BABAR detector at the PEP-II storage ring during 2000 and correspond to an integrated luminosity of 20.3 $`\mathrm{fb}^1`$ collected at the $`\mathrm{{\rm Y}}(4S)`$ resonance and 2.6 $`\mathrm{fb}^1`$ collected 40 MeV below the resonance for background studies (off-peak events).

PEP-II is an energy asymmetric storage ring, with positron and electron beam energies of about 3.11 and 9.0 GeV. The center-of-mass frame of the $`e^+e^{}`$ collision is therefore boosted along the $`z`$ direction in the laboratory frame, enabling decay time-dependent measurements of $`B`$ mesons through vertex reconstruction.

Samples of simulated $`B\overline{B}`$ and continuum events were analysed through the same analysis chain as the data. The equivalent luminosity of the generic simulated data is approximately equal to one third of the on-resonance data, while the equivalent luminosity of a specialized simulated sample containing $`B^0D^{}\pi ^+`$ decays followed by $`D^{}\overline{D}{}_{}{}^{0}\pi _{}^{}`$ is 2.9 times larger than the on-resonance data sample.

A detailed description of the BABAR detector and the algorithms used for the track reconstruction, particle identification and selection of $`B\overline{B}`$ events is provided elsewhere ; a brief summary is given here.

Only charged particles are used for the partial reconstruction of the signal. Particles with momentum higher than $`170\mathrm{MeV}/c`$ are reconstructed by matching hits in the silicon vertex tracker (SVT) with track elements in the drift chamber (DCH). Since tracks with momentum below $`170\mathrm{MeV}/c`$ do not leave signals on many wires in the DCH due to the bending induced by the magnetic field, they are reconstructed in the SVT alone.

Electron and muon identification is used in veto mode for the selection of the fast pion. Electrons are identified on the basis of the energy deposited in the electromagnetic calorimeter (EMC), the track momentum and the energy loss in the DCH. Muons are selected by requiring deep penetration in the instrumented flux return (IFR).

The Cherenkov light emission measured in the particle identification detector (DIRC) is employed to reject kaons from the fast pion sample.

Neutral particles are reconstructed from clusters in the EMC that are unmatched to projected charged tracks. They are used only to compute event shape quantities.

## 3 The partial reconstruction technique

The $`B^0D^{}\pi ^+`$ decays are reconstructed using only the tracks from the $`\pi ^+`$ (fast pion) and from the $`\pi ^{}`$ (slow pion) in $`D^{}\overline{D}{}_{}{}^{0}\pi _{}^{}`$. The fast pion momentum in the $`\mathrm{{\rm Y}}(4S)`$ rest frame is required to be in the range kinematically allowed for the decay $`B^0D^{}\pi ^+`$ . For nominal values of the beam energies this momentum is in the range 2.114 - 2.404 GeV/$`c`$. The momentum of the slow pion is required to be greater than 50 MeV/$`c`$.

Assuming that the fast and slow pion come from the decay $`B^0D^{}\pi ^+`$ followed by the two body decay $`D^{}X\pi ^{}`$, it is possible to compute the mass of the recoiling system $`X`$ averaged over an unmeasured angle due to the unknown direction of the $`B`$ momentum in the $`\mathrm{{\rm Y}}(4S)`$ rest frame. For $`B^0D^{}\pi ^+`$ decays this recoil mass $`M_{rec}`$ peaks at the $`D^0`$ mass with a width slightly less than 3 $`\mathrm{MeV}/c^2`$. Assuming that the $`B`$ momentum lies in the plane defined by the fast and slow pions momenta in the $`\mathrm{{\rm Y}}(4S)`$ rest frame, it is possible to compute the helicity angle of the pion in the $`D^{}`$ rest frame as well as the candidate $`\overline{D}^0`$ direction.

Since the dominant source of background is continuum events, the selection procedure aims at reducing this contribution. The main requirements of this selection are the following:

* R2, the ratio of the second to the zeroth Fox-Wolfram moment , computed from charged particles, is required to be less than 0.35.
* No other tracks should be in a cone of opening angle 0.4 rad centered on the fast pion momentum in the $`\mathrm{{\rm Y}}(4S)`$ rest frame. This cut is effective against continuum events because in this case tracks tend to be clustered in jets.
* A Fisher discriminant $`F_D`$ is computed from 15 event shape variables. Among these are the scalar sum of the momenta in the $`\mathrm{{\rm Y}}(4S)`$ rest frame of all tracks and neutrals, in nine $`20^0`$ angular bins around the fast pion direction. The signal peaks at $`F_D=0.5`$ while the continuum background peaks at $`F_D=0`$. Events are required to satisfy $`F_D<0.1`$.
* The cosine of the soft pion helicity angle is required to be larger than 0.4 in absolute value.

The $`B^0`$ decay point is determined from a vertex fit of the fast and slow pion tracks and the beam spot position in the plane perpendicular to the beam axis (the $`x`$-$`y`$ plane). The beam spot is determined on a run-by-run basis using two-prong events . Its size in the horizontal direction is 120 $`\mu `$m. Although the beam spot size in the $`y`$ direction is only a few $`\mu `$m, a beam spot constraint of 30 $`\mu `$m is applied to account for the flight of the $`B^0`$ in the $`y`$ direction. Only events for which the probability of the vertex fit is greater than 0.1% are considered further.

The decay point of the other $`B`$ is determined from a subset of the remaining tracks in the event. All the tracks with a center-of-mass angle greater than 1 rad with respect to the candidate $`\overline{D}^0`$ direction are considered. This requirement is used to remove most of the tracks from the decay of the $`\overline{D}^0`$ daughter of the $`D^{}`$, which would otherwise bias the reconstruction of the other /B/ vertex position. The selected tracks are then constrained to the beam spot in the $`x`$-$`y`$ plane. The track with the largest contribution to the vertex $`\chi ^2`$, if greater than 6, is removed and the fit iterated until no tracks fails this requirement. Vertices composed of just one track and the beam spot are rejected in order to reduce the number of poorly measured vertices. Simulation shows that after all these requirements in about 85% of signal events the other vertex has no tracks from the $`\overline{D}^0`$ decay.

Figure 1 shows the recoil mass distribution for on-resonance data events obtained by applying all the above cuts when the fast and slow pions have opposite-sign charges. These events are used to measure the $`B^0`$ lifetime. Events with same-sign charges for the two pions are used as a background control sample (same-sign sample). Another background control sample is obtained by reversing the cut on the Fisher discriminant and requiring $`F_D>0`$ ($`B\overline{B}`$ depleted sample).

The individual distributions shown in Figure 1 are obtained by fitting simultaneously to the on- and off-resonance data the contributions from

* $`B^0D^{}\pi ^+`$ events;
* $`B^0D^{}\rho ^+`$ events;
* $`B^0D^{}\mathrm{}^+\nu `$ events;
* peaking $`B\overline{B}`$ background, where the fast and slow pions originate from the same B meson decay;
* non-peaking $`B\overline{B}`$ background, excluding the events from the previous four categories;
* peaking $`c\overline{c}`$ background, where the slow pion comes from the $`D^{}\overline{D}{}_{}{}^{0}\pi _{}^{}`$ decay;
* continuum background excluding the events from the previous category.

All peaking backgrounds peak in the recoil mass distribution.

The shape of each contribution has been determined from the Monte Carlo simulation while all the normalizations have been left free in the fit except those for $`B^0D^{}\rho ^+`$ and $`B^0D^{}\mathrm{}^+\nu `$ events which have been fixed to the value from the Monte Carlo simulation. The $`B^0D^{}\rho ^+`$ branching fraction is taken to be $`(6.8\pm 3.4)10^3`$  giving an estimate of 1570 events in the final sample from this decay mode. The fit gives $`6971\pm 241`$ $`B^0D^{}\pi ^+`$ events. In the lifetime analysis to be described later, the $`B^0D^{}\rho ^+`$ events will be included as signal in addition to the 6971 $`B^0D^{}\pi ^+`$ events. Table 1 reports the composition of the selected sample for recoil mass larger than 1.86 $`GeV/c^2`$ (signal region sample) as determined by this fit.

## 4 The lifetime measurement

The PEP-II collider produces $`B\overline{B}`$ pairs moving along the beam direction ($`z`$ axis) with an average Lorentz boost of $`\beta \gamma =0.55`$. The lifetime is determined by measuring the quantity $`\mathrm{\Delta }z=z_{decay}z_{other}`$, where $`z_{decay}`$ ($`z_{other}`$) is the position along the beam line of the reconstructed $`B^0D^{}\pi ^+`$ decay (other) vertex. To remove badly reconstructed vertices, all events for which $`\sigma _{\mathrm{\Delta }z}>400\mu \mathrm{m}`$, where $`\sigma _{\mathrm{\Delta }z}`$ is the uncertainty on $`\mathrm{\Delta }z`$ computed for each event, are rejected.

Residual tracks from the $`\overline{D}^0`$ decay, not removed by the track selection for the other $`B`$ vertex, bias the reconstruction of the other $`B`$ vertex position. This bias is removed by correcting the $`\mathrm{\Delta }z`$ value for each event with a correction function determined from the simulated signal sample as a function of $`\mathrm{\Delta }z`$.

The proper time difference between $`B`$ decays is then computed with the relation $`\mathrm{\Delta }t=\mathrm{\Delta }z/c\beta \gamma `$. A fit with a double Gaussian to the $`\mathrm{\Delta }t`$ residuals in the Monte Carlo simulation shows that 75% of the events are contained in the narrower Gaussian, which has a width of 0.8 $`\mathrm{ps}`$.

The lifetime $`\tau _{B^0}`$ is obtained from an unbinned maximum likelihood fit to the two-dimensional $`\mathrm{\Delta }t`$, $`\sigma _{\mathrm{\Delta }t}`$ distribution. The $`\mathrm{\Delta }t`$ distribution of signal events is described by the convolution of the decay probability distribution

$$f(\mathrm{\Delta }t_{true}|\tau _{B^0})=\frac{1}{2\tau _{B^0}}e^{|\mathrm{\Delta }t_{true}|/\tau _{B^0}},$$
(1)
with the experimental resolution function, which is represented by the sum of three Gaussian distributions. The first two, accounting for more than 99% of the events, have the form

$$G(\delta (\mathrm{\Delta }t),\sigma _{\mathrm{\Delta }t})=\frac{1}{\sqrt{2\pi }S\sigma _{\mathrm{\Delta }t}}e^{\frac{(\delta (\mathrm{\Delta }t)b\sigma _{\mathrm{\Delta }t})^2}{2S^2\sigma _{\mathrm{\Delta }t}^2}}$$
(2)
where $`\delta (\mathrm{\Delta }t)=\mathrm{\Delta }t\mathrm{\Delta }t_{true}`$ is the difference between the measured and the true value of $`\mathrm{\Delta }t`$, $`b`$ is a bias due to the charm tracks in the other vertex and the scale factor $`S`$ is introduced to account for possible misestimation of the error on the proper time difference $`\mathrm{\Delta }t`$. The third Gaussian of fixed bias $`b=0`$ and scale factor $`S=6`$ is added to account for badly mismeasured events (“outliers”).

The $`B^0D^{}\rho ^+`$ resolution function has been found with Monte Carlo simulation to be the same as for signal events. In the fit for the $`B^0`$ lifetime the $`B^0D^{}\rho ^+`$ events are considered as signal events. Events coming from the decay $`B^0D^{}\mathrm{}^+\nu `$ are added to the $`B\overline{B}`$ peaking component and the $`\mathrm{\Delta }t`$ distribution of the latter is assumed to be equal to the $`B\overline{B}`$ non-peaking component.

The $`\mathrm{\Delta }t`$ distribution of the B background contribution is described by a resolution function similar to that used for the signal. In this case the lifetime parameter represents an effective B lifetime which has been fitted independently from the signal lifetime.

The $`\mathrm{\Delta }t`$ distribution for $`uds`$ events is described by the sum of three Gaussians while the $`\mathrm{\Delta }t`$ distribution for $`c\overline{c}`$ events has been described by the sum of a Gaussian plus the same Gaussian convoluted with an exponential term to account for the effective charm hadron lifetime.

The function used to fit the data is the weighted sum of four contributions:

$`F(\mathrm{\Delta }t,\sigma ,\tau _{B^0})=[1f_{B\overline{B}}(M_{rec})f_{uds}(M_{rec})f_{c\overline{c}}(M_{rec})]F_{B^0}(\mathrm{\Delta }t,\sigma ,\tau _{B^0})+`$
$`f_{B\overline{B}}(M_{rec})F_{B\overline{B}}(\mathrm{\Delta }t,\sigma )+f_{uds}(M_{rec})F_{uds}(\mathrm{\Delta }t,\sigma )+f_{c\overline{c}}F_{c\overline{c}}(\mathrm{\Delta }t,\sigma )`$
where the functions $`F_{B^0}`$, $`F_{B\overline{B}}`$, $`F_{uds}`$ and $`F_{c\overline{c}}`$ describe the measured decay time difference distributions for the signal, $`B\overline{B}`$ background, $`uds`$ and $`c\overline{c}`$ events, respectively. $`f_{B\overline{B}}`$, $`f_{uds}`$ and $`f_{c\overline{c}}`$ are the probabilities that the event is from the $`B\overline{B}`$, $`uds`$ and $`c\overline{c}`$ background, computed for each event on the basis of the measured value of the recoil mass $`M_{rec}`$.

The key parameters describing the $`\mathrm{\Delta }t`$ distributions for the background events are fitted on the $`B\overline{B}`$ depleted control sample for the continuum background and on the same-sign sample for the $`B\overline{B}`$ background. The $`B\overline{B}`$ depleted control sample has been used because the integrated luminosity of the off resonance sample in not sufficient to precisely determine the parameters of the continuum background. An alternative parameterization obtained by fitting the side-band control sample has also been tried giving compatible results.

All the parameters describing the signal resolution function are free in the fit except the bias of the second Gaussian which is fixed to the value found on the signal Monte Carlo sample.

The result of the fit to the signal region sample in the range $`|\mathrm{\Delta }t|<15\mathrm{ps}`$ is $`\tau _{B^0}^{raw}=1.524\pm 0.040\mathrm{ps}`$, where the error is statistical only. Figure 2 shows the comparison between the measured $`\mathrm{\Delta }t`$ distribution and the fit result. All the fitted parameters appear in Table 2.

The raw value $`\tau _{B^0}^{raw}`$ value is corrected for a small bias, $`0.014\pm 0.020\mathrm{ps}`$, observed when fitting the signal Monte Carlo sample using the same procedure, yielding the corrected value $`\tau _{B^0}=1.510\pm 0.040\mathrm{ps}`$.

## 5 Systematic uncertainties and cross-checks

The systematic error on $`\tau _{B^0}`$ is computed by adding in quadrature the contribution from several sources, described below and summarised in Table 3.

The fractions of the various background components are varied by their uncertainties obtained from the recoil mass fit. The parameters of the various background $`\mathrm{\Delta }t`$ distributions are also varied by their uncertainties, properly accounting for their mutual correlations. The effective lifetime of the $`B\overline{B}`$ peaking component is varied by $`\pm 0.044\mathrm{ps}`$ corresponding to the difference between the value found in the data and in the Monte Carlo sample. An alternate resolution function for this component has also been tried with negligible variation on $`\tau _{B^0}`$.

The bias of the wide Gaussian, the only parameter of the signal resolution function which is not fitted, is varied in a conservative range. Several different analytical expressions are used to represent the small fraction of outliers. The fit range is varied from ($`10`$, $`10`$) $`\mathrm{ps}`$ to ($`20`$, $`20`$) $`\mathrm{ps}`$. The parameters of the $`\mathrm{\Delta }z`$ correction are varied according to the uncertainty due to the finite signal Monte Carlo sample size. The systematic uncertainty related to this correction is estimated from the effect of a $`\pm 5\%`$ variation on the fraction of fitted other-$`B`$ vertices which have no tracks from the $`\overline{D}^0`$ decay.

The $`z`$ length scale is determined to about $`0.4\%`$ from secondary interactions with a beam pipe section of known length. The statistical uncertainty of the residual bias found on signal Monte Carlo events is also added to the systematic error.

The total systematic error of $`\pm 0.038\mathrm{ps}`$ is found by adding in quadrature the uncertainties from the above sources.

The dependence of the result on several different variables (angular width of the cone used to reject the $`\overline{D}^0`$ tracks, fast pion momentum, polar and azimuthal angles) has been studied: no statistically significant effect is found.

## 6 Conclusion

In conclusion the neutral $`B`$ meson lifetime has been measured with a sample of $`6971\pm 241`$ partially reconstructed $`B^0D^{}\pi ^+`$ decays:

$$\tau _0=1.510\pm 0.040\pm 0.038\mathrm{ps}.$$

This preliminary value is consistent with other recent BABAR measurements  and with the world average $`B^0`$ lifetime .

## 7 Acknowledgments

We are grateful for the extraordinary contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible. The success of this project also relies critically on the expertise and dedication of the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and the kind hospitality extended to them. This work is supported by the US Department of Energy and National Science Foundation, the Natural Sciences and Engineering Research Council (Canada), Institute of High Energy Physics (China), the Commissariat à l’Energie Atomique and Institut National de Physique Nucléaire et de Physique des Particules (France), the Bundesministerium für Bildung und Forschung (Germany), the Istituto Nazionale di Fisica Nucleare (Italy), the Research Council of Norway, the Ministry of Science and Technology of the Russian Federation, and the Particle Physics and Astronomy Research Council (United Kingdom). Individuals have received support from the A. P. Sloan Foundation, the Research Corporation, and the Alexander von Humboldt Foundation.


text:
# Exact results for a tunnel-coupled pair of trapped Bose-Einstein condensates

## Abstract

A model describing coherent quantum tunneling between two trapped Bose-Einstein condensates is shown to admit an exact solution. The spectrum is obtained by the algebraic Bethe ansatz. An asymptotic analysis of the Bethe ansatz equations leads us to explicit expressions for the energies of the ground and first excited states in the limit of weak tunneling and all energies for strong tunneling. The results are used to extract the asymptotic limits of the quantum fluctuations of the boson number difference between the two Bose-Einstein condensates and to characterize the degree of coherence in the system.

Macroscopic quantum tunneling is one of the most fascinating phenomena in condensed matter physics . Experimental realization of Bose-Einstein condensates (BEC’s) in dilute atomic alkali gases has stimulated a diverse range of theoretical and experimental research activity . A particularly exciting possibility is that a pair of BEC’s (such as a BEC trapped in a double-well potential) may provide a model tunable system in which to observe macroscopic quantum tunneling. In this Letter we show that the canonical Hamiltonian for a pair of tunnel-coupled BEC’s  has an exact solution. The model is also realizable in Josephson-coupled super-conducting metallic nano-particles . This connects the model to the rich and powerful mathematics associated with the algebraic Bethe ansatz  and provides a means to study the degradation of the coherence between the two BEC’s which occurs as the tunnel coupling increases (or the size of the BEC decreases). The exactness of our approach means that it is not necessary to resort to approximations such as Gross-Pitaevskii mean-field theory and the phase-number formulation which become questionable in the coherent-incoherent crossover region.

The canonical Hamiltonian takes the form

$`H`$ $`=`$ $`{\displaystyle \frac{K}{8}}(N_1N_2)^2{\displaystyle \frac{\mathrm{\Delta }\mu }{2}}(N_1N_2)`$ (2)
$`{\displaystyle \frac{_𝒥}{2}}(a_1^{}a_2+a_2^{}a_1).`$
where $`a_1^{},a_2^{}`$ denote the single-particle creation operators in the two wells and $`N_1=a_1^{}a_1,N_2=a_2^{}a_2`$ are the corresponding boson number operators. The total boson number $`N_1+N_2`$ is conserved and set to the fixed value of $`N`$. The physical meaning of the coupling parameters for different realizable systems may be found in Ref. . It is useful to divide the parameter space into three regimes; viz. Rabi ($`K/_J<<N^1`$), Josephson ($`N^1<<K/_J<<N`$) and Fock ($`N<<K/_J`$). In the Josephson region one expects coherent superposition of the two BEC’s (Schrödinger cat states) to be possible whereas in the Fock region the two BEC’s will be (in some sense) localised in the two separate wells. There is a correspondence between (2) and the motion of a pendulum . In the Rabi and Josephson regimes this motion is semiclassical, unlike the case of the Fock regime. For both the Fock and Josephson regimes the analogy corresponds to a pendulum with fixed length, while in the Rabi regime the length varies. An important problem is to study the behaviour in the crossover regimes.

The exact solvability of (2) which we discuss here follows from the fact that it is mathematically equivalent to the discrete self-trapping dimer model studied by Enol’skii et al. , who solved the model through the algebraic Bethe ansatz. The Lax operator appearing in the Yang-Baxter equation which underlies the exact solution first appeared in . The consequences of this result are legion. The exact solution gives a direct method, through numerical analysis, to investigate the nature of the system in the crossover regions between the Fock, Rabi and Josephson regimes. It also provides a means to test the limits of applicability of previous approximate treatments . Solvability through the Yang-Baxter equation raises questions about the probability distribution of the energy level spacings . Moreover, the formulation of the model through the algebraic version of the Bethe ansatz opens possibilities for the calculation of time-dependent form factors and correlation functions, as achieved in  for $`\mathrm{\Delta }\mu =0`$.

In this Letter we will show that the exact solution allows for the analysis of the asymptotic behaviour of the energy spectrum in the limits of strong and weak tunneling. These results are used to extract the asymptotic limits of the quantum fluctuations of the relative particle number between the two BEC’s and the degree of coherence in the system for the ground state. Our approach also makes it feasible to undertake an asymptotic analysis of the system at finite temperature which gives insight into the relative influence of thermal and quantum fluctuations (cf. ). We mention that these findings provide a useful tool in the numerical evaluation of the exact solution in the strong and weak tunneling regions. In contrast to most Bethe ansatz solvable models, the model discused here, as in the case of the Azbel-Hofstadter problem, is a quantum mechanical model rather than a quantum field theory; i.e., it has a fixed and finite number of degrees of freedom.

Bethe ansatz solution. Following the standard procedure of the algebraic Bethe ansatz , we can derive the Bethe ansatz equations (BAE). For $`N`$ total bosons the BAE for (1) can be written in the form

$$\eta ^2(v_\alpha ^2\omega ^2)=\underset{\beta \alpha }{\overset{N}{}}\frac{v_\alpha v_\beta \eta }{v_\alpha v_\beta +\eta },\alpha =1,\mathrm{},N.$$

The parameters $`\eta ,\omega `$ which naturally arise in the Bethe ansatz solution, along with a scaling factor $`\kappa `$, are related to the coupling constants of (2) through the identification

$$K=2\kappa \eta ^2,\mathrm{\Delta }\mu =2\kappa \eta \omega ,_𝒥=2\kappa .$$

Each set of numbers $`\{v_\alpha \}_{\alpha =1}^N`$ which is a solution of the BAE defines an eigenstate of the Hamiltonian. For details of the derivation of (2) we refer to . For such a solution, the eigenvector has the form

$$\mathrm{\Psi }(v_1,\mathrm{},v_N)=\underset{\alpha =1}{\overset{N}{}}C(v_\alpha )|0$$

with $`C(u)=(u\omega +\eta N_2)a_1^{}+\eta ^1a_2^{}`$ and $`|0`$ the Fock vacuum.

For each $`N`$ we expect $`N+1`$ independent solutions of the BAE. Note that in the derivation of the BAE it is assumed that $`v_\alpha `$ are distinct for different $`\alpha `$. This is a result of the Pauli principle for Bethe ansatz solvable models as proved by Izergin and Korepin  for the one-dimensional $`\delta `$-function interaction boson gas. Their result may be adapted to the present case, which plays an important role in our asymptotic analysis of the Bethe ansatz solutions. From the BAE, we may derive the useful identity

$`{\displaystyle \underset{\alpha =1}{\overset{m}{}}}\eta ^2(v_\alpha \omega ^2)={\displaystyle \underset{\alpha =1}{\overset{m}{}}}{\displaystyle \underset{\beta =m+1}{\overset{N}{}}}{\displaystyle \frac{v_\alpha v_\beta \eta }{v_\alpha v_\beta +\eta }}`$ (3)
which will be used frequently.

For a given solution to the Bethe ansatz equations, the corresponding energy eigenvalue of the Hamiltonian is

$`E`$ $`=`$ $`\kappa (\eta ^2{\displaystyle \underset{\alpha }{\overset{N}{}}}(1+{\displaystyle \frac{\eta }{v_\alpha u}}){\displaystyle \frac{\eta ^2N^2}{4}}u\eta Nu^2`$ (5)
$`\eta ^2+\omega ^2+(u^2\omega ^2){\displaystyle \underset{\alpha }{\overset{N}{}}}(1{\displaystyle \frac{\eta }{v_\alpha u}})).`$
Note that this expression is independent of the spectral parameter $`u`$ which can be chosen arbitrarily. The formula simplifies considerably with the choice $`u=\omega `$, by employing (3), which yields a polynomial form. However, for the purpose of an asymptotic analysis in the Rabi regime, it is more convenient to choose $`u=0`$, while for the Fock regime we use $`u=\eta ^2`$.

Asymptotics. We start our analysis with the Rabi regime where $`\eta ^2N<<1`$. From the BAE it is clear that $`\eta ^2v_\alpha ^21`$ as $`\eta 0`$, so that $`v_\alpha \pm \eta ^1`$. However, when $`\eta =0`$ we know that the Hamiltonian is diagonalizable by using the Bogoliubov transformation, from which we can deduce that the solution of the BAE corresponding to the ground state must have $`v_\alpha \eta ^1`$. Excitations correspond to changing the signs of the leading terms in the Bethe ansatz roots. To study the asymptotic behaviour for the $`m`$th excited state, we set

$`v_\alpha `$ $``$ $`\eta ^1+ϵ_\alpha +\eta \delta _\alpha ,\alpha =1,\mathrm{},m,`$ (6)
$`v_\alpha `$ $``$ $`\eta ^1+ϵ_\alpha +\eta \delta _\alpha ,\alpha =m+1,\mathrm{},N,`$ (7)
with the convention that the ground state corresponds to $`m=0`$.

From the leading terms of the BAE for $`v_\alpha ,\alpha m`$ we find

$$ϵ_\alpha =\underset{\beta \alpha }{\overset{m}{}}\frac{1}{ϵ_\alpha ϵ_\beta },$$
(8)
which implies

$$\underset{\alpha =1}{\overset{m}{}}ϵ_\alpha =0,\underset{\alpha =1}{\overset{m}{}}ϵ_\alpha ^2=\frac{m(m1)}{2}.$$

In a similar fashion we have for $`m<\alpha N`$

$$ϵ_\alpha =\underset{\alpha \beta =m+1}{\overset{N}{}}\frac{1}{ϵ_\alpha ϵ_\beta },$$
(9)
which implies

$$\underset{\alpha =m+1}{\overset{N}{}}ϵ_\alpha =0,\underset{\alpha =m+1}{\overset{N}{}}ϵ_\alpha ^2=\frac{(Nm)(Nm1)}{2}.$$

It is clear from (8) and (9) why the Pauli exclusion principle applies in the present case. In the asymptotic expansion for $`v_\alpha `$, $`ϵ_\alpha `$ is assumed finite. However, if $`v_\alpha =v_\beta `$ for some $`\alpha ,\beta `$, then $`ϵ_\alpha =ϵ_\beta `$ and (8) and (9) imply that $`ϵ_\alpha ,ϵ_\beta `$ are infinite which is a contradiction. Hence $`v_\alpha `$ must be distinct for different $`\alpha `$. Note also that for this approximation to be valid we require $`\eta ^1>>ϵ_\alpha .`$ However, we see that $`|ϵ_\alpha |`$ is of the order of $`N^{1/2}`$. Thus our approximation will be valid for $`\eta N^{1/2}<<1`$, which is precisely the criterion for the Rabi region and consequently $`N`$ cannot be arbitrarily large for fixed $`\eta `$, or vice versa.

Now we go to the next order. From (3) we find

$`{\displaystyle \underset{\alpha =1}{\overset{m}{}}}\delta _\alpha `$ $`=`$ $`{\displaystyle \frac{m(m1)}{4}}+{\displaystyle \frac{m(mN)}{2}}{\displaystyle \frac{m\omega ^2}{2}}`$ (10)
$`{\displaystyle \underset{\alpha =m+1}{\overset{N}{}}}\delta _\alpha `$ $`=`$ $`{\displaystyle \frac{(Nm)(Nm1)}{4}}+{\displaystyle \frac{m(mN)}{2}}`$ (12)
$`+{\displaystyle \frac{(Nm)\omega ^2}{2}}`$
which using (5) leads us to the result

$$\frac{E_m}{\kappa }N+2m\frac{\eta ^2\omega ^2(N2m)}{2}+\frac{\eta ^2N}{4}+\frac{\eta ^2}{2}m(Nm).$$

The energy level spacings $`\mathrm{\Delta }_m=E_mE_{m1}`$ are thus

$`\mathrm{\Delta }_m`$ $``$ $`\kappa \left(2+\eta ^2\omega ^2+{\displaystyle \frac{\eta ^2}{2}}(N2m+1)\right).`$ (13)
One may check that $`\mathrm{\Delta }_m/N`$ is of the order of $`N^1`$. This indicates that the Rabi regime is semiclassical. Our value for the gap between the ground and first excited state agrees, to leading order in $`\eta ^2N`$, with the Gross-Pitaevskii mean-field theory which gives a Josephson plasma frequency of $`\omega _J=2\kappa (1+\eta ^2N/2)^{1/2}`$.

Now we look at the asymptotic behaviour of the BAE in the Fock regime $`\eta ^2>>N`$. It is necessary to distinguish the following cases: (i) $`\omega =0`$ and (ii) $`\omega 0`$.

(i) $`\omega =0`$. In this case, it is appropriate to consider the permutation operator $`P`$ which interchanges the labels 1 and 2 in (2). For $`\omega =0`$, $`P`$ commutes with the Hamiltonian, and any eigenvector of the Hamiltonian is also an eigenvector of $`P`$ with eigenvalue $`\pm 1`$. Therefore the Hilbert space splits into the direct sum of two subspaces corresponding to the symmetric and antisymmetric eigenfunctions. From now on we restrict ourselves to the case when $`N`$ is even, i.e., $`N=2M`$, although a similar calculation is also applicable to the case when $`N`$ is odd. A careful analysis leads us to conclude that the ground state lies in the symmetric subspace. The asymptotic form of the roots of the BAEs for the ground state takes the “string”-like structure

$`v_{j\pm }`$ $``$ $`(Mj)\eta \pm i{\displaystyle \frac{C_M^j}{(j1)!}}\eta ^{(2j1)}`$ (15)
$`+M(M+1)\eta ^3\delta _{j1},j=1,\mathrm{},M.`$
where $`C_M^j`$ is a binomial coefficient. For this asymptotic ansatz to be valid, we require that any term in the asymptotic expansion should be much smaller than those preceeding. This yields $`\eta ^2>>N`$ which coincides with the defining condition for the Fock region. Throughout, the Pauli exclusion principle has been taken into account to exclude any possible spurious solutions of the BAE.

The above structure clearly indicates that in the ground state the $`N`$ bosons fuse into $`M`$ “bound” states and excitations correspond to a breakdown of these bound states. Specifically, the first and second excited states correspond to the breakdown of the bound state at $`(M1)\eta `$, with the first excited state in the antisymmetric subspace and the second excited state in the symmetric subspace. Explicitly, we can write down the spectral parameter configurations for the first two excited states

$`v_{1+}`$ $``$ $`M\eta +a_{1+}\eta ^3,v_1(M1)\eta +a_1\eta ^3,`$ (16)
$`v_{j\pm }`$ $``$ $`(Mj)\eta +a_{j\pm }\eta ^{(2j1)},j=2,\mathrm{},M,`$ (17)
with

$`a_{1+}`$ $``$ $`{\displaystyle \frac{M+1}{2}},a_1{\displaystyle \frac{M(M+1)}{2}},`$ (18)
$`a_{2\pm }`$ $``$ $`{\displaystyle \frac{(M1)^2\pm (M1)\sqrt{13M^2+10M+1}}{12}},`$ (19)
$`a_{3\pm }`$ $``$ $`\pm {\displaystyle \frac{(M1)(M2)\sqrt{2M(M+1)}}{24}},`$ (20)
$`a_{j\pm }`$ $``$ $`{\displaystyle \frac{Mj+1}{\sqrt{(j+1)j(j1)(j2)}}}a_{j1,\pm },j=3,\mathrm{},M,`$ (21)
for the (antisymmetric) first excited state and

$`a_{1+}`$ $``$ $`{\displaystyle \frac{(M+1)(2M+1)}{2}},a_1{\displaystyle \frac{M(M+1)}{2}},`$ (22)
$`a_{2\pm }`$ $``$ $`{\displaystyle \frac{(M1)^2\pm i(M1)\sqrt{11M^2+14M1}}{12}},`$ (23)
$`a_{3\pm }`$ $``$ $`\pm i{\displaystyle \frac{(M1)(M2)\sqrt{2M(M+1)}}{24}},`$ (24)
$`a_{j\pm }`$ $``$ $`{\displaystyle \frac{Mj+1}{\sqrt{(j+1)j(j1)(j2)}}}a_{j1,\pm },j=3,\mathrm{},M,`$ (25)
for the (symmetric) second excited state. The breakdown of the bound state at $`(Mj)\eta `$, $`j=2,\mathrm{},M`$ results in the higher excited states.

Substituting these results into (5) leads us to the asymptotic ground state energy

$$E_02\kappa \eta ^2M(M+1),$$

while for the first and second excited states we have

$`E_1`$ $``$ $`\kappa \eta ^2\kappa \eta ^2{\displaystyle \frac{M^2+M2}{3}},`$ (26)
$`E_2`$ $``$ $`\kappa \eta ^2+\kappa \eta ^2{\displaystyle \frac{5M^2+5M+2}{3}}.`$ (27)
In contrast to the Rabi regime, the Fock regime is not semiclassical, as the ratio of the gap $`\mathrm{\Delta }`$ and $`N`$ is of finite order when $`N`$ is large.

We can perform a similar analysis for odd $`N`$. In this case, the gap between the ground and the first excited states is proportional to $`\kappa \eta ^2`$ instead of $`\kappa \eta ^2`$. This indicates there is a strong parity effect in the Fock regime. Its possible physical implication remains to be explored.

(ii) $`\omega 0`$. In this case the root structure is somewhat more complicated than for $`\omega =0`$, so we will not present the details. We remark however that our calculations show that up to order $`\eta ^2`$ the ground state energy eigenvalue takes the same form as in the case $`\omega =0`$. Actually, the leading contribution arising from the $`\omega `$ term appears only as $`\omega ^2\eta ^4`$. This means that the results presented below are applicable for all values of $`\omega `$ (or equivalently $`\mathrm{\Delta }\mu `$).

Quantum fluctuations and coherence factor. Although it is difficult to define rigorously, the relative phase between Bose-Einstein condensates has been useful in understanding interference experiments. By definition, the relative phase $`\mathrm{\Phi }`$ is conjugate to the relative number of atoms in the two condensates $`nN_1N_2`$. Using the Feynmann-Hellman theorem, we find that

$$<\mathrm{\Delta }n^2>=8\frac{E}{K}4(\frac{E}{\mathrm{\Delta }\mu })^2.$$

For the ground state in the limit of strong  tunneling (i.e., Rabi regime), $`<\mathrm{\Delta }n^2>=N(\mathrm{\Delta }\mu N/_J)^2`$. In the case of weak  tunneling (i.e., Fock regime), $`<\mathrm{\Delta }n^2>2N(N+2)(_J/K)^2`$. The degree of coherence between the two BEC’s can be discussed in terms of

$$\alpha \frac{1}{2N}<a_1^{}a_2+a_2^{}a_1>=\frac{1}{N}\frac{E}{_J}$$

where the last identity follows from the Feynmann-Hellmann theorem. In terms of the relative phase $`\mathrm{\Phi }`$, $`\alpha <\mathrm{cos}\mathrm{\Phi }>`$ . In the strong coupling limit, $`\alpha 1N^1(\mathrm{\Delta }\mu )^2/(8_J)^2`$, indicating very close to full coherence in the ground state. In the opposite limit, we have $`\alpha 2(N+2)_J/K<<1`$, indicating the absence of coherence.

The above results give the first order corrections to the results presented in  for the number fluctuations and the coherence factor at zero temperature. Moreover, we can extend this approach to gain insight into the effects of thermal fluctuations. This problem was recently studied in  using Gross-Pitaevskii theory. In the Fock regime, the energy gap is large, and so for the low temperature thermodynamics it is sufficient to consider the first two excited states, for which the calculations are straightforward. However, in the Rabi regime the small energy scale requires us to consider the entire spectrum. The asymptotic form of the partition function is found to be (valid for all temperatures)

$`Z`$ $`=`$ $`{\displaystyle \underset{m}{}}\mathrm{exp}(E_m/T)`$ (28)
$``$ $`{\displaystyle \frac{\mathrm{sinh}\frac{\beta }{2}(N+1)}{\mathrm{sinh}\frac{\beta }{2}}}+{\displaystyle \frac{\beta K}{4_J}}\mathrm{exp}(\beta N)\widehat{𝒟}\left({\displaystyle \frac{1\mathrm{exp}\beta (N+1)}{1\mathrm{exp}\beta }}\right),`$ (29)
with $`\beta =\frac{2\kappa }{T}`$ (we set Boltzmann’s constant $`k_B=1`$) and $`\widehat{𝒟}`$ is the differential operator

$`\widehat{𝒟}`$ $``$ $`{\displaystyle \frac{(\mathrm{\Delta }\mu )^2N}{4_JK}}+{\displaystyle \frac{N}{4}}+\left({\displaystyle \frac{(\mathrm{\Delta }\mu )^2}{2_JK}}+{\displaystyle \frac{N}{2}}\right){\displaystyle \frac{}{\beta }}{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{^2\beta }}.`$ (30)
In the usual way it is possible to derive a number of physical quantities from $`Z`$. For example, the temperature dependent coherence factor is simply

$$\alpha (T)=\frac{T}{N}\frac{\mathrm{ln}Z}{_J}$$

and an expression for the particle number fluctuations can be found in a similar manner. This provides an analytic means to investigate the thermal effects for the Rabi regime. In order to investigate the crossover to weaker tunneling, it is possible to solve the BAE numerically to compute these quantities. This will be the subject of future work.

In conclusion, we have shown that the canonical Hamiltonian describing quantum tunneling between a pair of BEC’s is integrable in the sense of the algebraic Bethe ansatz. An closed form expression was given for all of the energy eigenstates and eigenvalues in terms of the solutions of the Bethe ansatz equations. Our approach makes the calculation of form factors  and correlation functions possible without any of the limitations associated with Gross-Pitaevskii mean-field theory or the number-phase representation.

We thank J. Corney, G.J. Milburn, and M.D. Gould for helpful discussions. This work was supported by the Australian Research Council.


text:
# 1 Introduction

## 1 Introduction

The $`ep`$ collider HERA is a unique accelerator because it enables us to investigate the proton structure in extremely deep region. Furthermore we can study the basic electron-quark interactions at the high energies only by HERA which provides us both baryonic and leptonic quantum numbers in the initial states. In fact HERA could set the best limit on the masses and the coupling strengths of leptoquarks via detailed analyses of the deep inelastic (DIS) neutral current (NC) and charged current (CC) processes.

It is widely known that the scalar partners of quarks (squarks : $`\stackrel{~}{q}`$) have similar properties to the leptoquaks in a class of the supersymmetric (SUSY) standard models, the R-parity breaking (RB) models . Here the R-parity is defined by $`R=(1)^{L+3B+2S}`$, where $`L`$, $`B`$ and $`S`$ denote the lepton number, the baryon number and the spin, respectively. The models are characterized by the R-parity violating $`\widehat{L}\widehat{Q}\widehat{D}`$ type superpotentials, from which the leptoquark interactions $`\mathrm{}q\stackrel{~}{q}`$ are derived. The superfields $`\widehat{L}`$, $`\widehat{Q}`$ and $`\widehat{D}`$ contain, respectively, the left-handed leptons, the left-handed quarks and the right-handed down-quark together with the SUSY partners $`\stackrel{~}{\mathrm{}}_L`$, $`\stackrel{~}{q}_L`$ and $`\stackrel{~}{d}_R`$. We could study in detail the single production of squarks as well as leptoquarks by HERA experiments.

In 1997 both H1  and ZEUS  collaborations reported an event excess at large $`x`$ and high $`Q^2`$ in comparison with the Standard Model (SM) expectations in the DIS $`e^+pe^+X`$. The news very much excited the high energy physics community. Various ideas to understand the anomaly have extensively been examined by theoreticians since then . We have also proposed an interpretation of the anomalous events by the scalar top quark (stop) production in the SUSY models with RB interactions . Contrary to initial expectation, the anomalies have gradually faded away from the whole data sample with increasing experimental data. However, this fact impressed us that HERA could have the potentiality exploring physics beyond the SM.

In this paper we consider another type of the R-parity violating SUSY model, with the $`\widehat{L}\widehat{L}\widehat{E}`$ type superpotential, where $`\widehat{E}`$ contains the right-handed leptons together with the SUSY partners $`\stackrel{~}{\mathrm{}}_R`$. We show that the multilepton events come from the single production of a scalar partner of the neutrinos (sneutrinos) at HERA. As the process is not an $`eq`$ scattering of the resonant type, the number of events is too small to be detected at the present HERA with the integrated luminosity of around $`100`$pb-1. Fortunately, however, HERA has been upgraded, and experiments using polarized $`e^{}`$ or $`e^+`$ beams with high luminosity will soon start. Then it will be expected that we are able to search for rare events with higher statistics than before.

## 2 Models and constraints

In the minimal SUSY standard model (MSSM), the general RB superpotential $`W_{/R}`$ is written by

$$W_{/R}=\lambda _{ijk}\widehat{L}_i\widehat{L}_j\widehat{E^c}_k+\lambda _{ijk}^{}\widehat{L}_i\widehat{Q}_j\widehat{D^c}_k+\lambda _{ijk}^{\prime \prime }\widehat{U^c}_i\widehat{D^c}_j\widehat{D^c}_k,$$
(1)
where $`i`$, $`j`$ and $`k`$ are generation indices. The first two terms violate the lepton number $`L`$ and the last term violates the baryon number $`B`$. Incorporating RB interactions into the MSSM we have a possibility to unveil yet unresolved problems as (i) the cosmic baryon number violation, (ii) the origin of the masses and the magnetic moments of neutrinos and (iii) some interesting rare processes induced by the $`L`$ and/or $`B`$ violation. Here we consider the first $`\widehat{L}_i\widehat{L}_j\widehat{E}_k`$ term, in which we set $`(i,j,k)=(2,3,1)`$. The Yukawa-type interaction Lagrangian

$$L=\lambda _{231}(\overline{\stackrel{~}{e_R}}\overline{\nu _\mu }\tau _L\overline{\stackrel{~}{e_R}}\overline{\nu _\tau }\mu _L\stackrel{~}{\nu _\tau }\overline{e^c}\mu _L+\mathrm{})+h.c.$$
(2)
is derived from the superpotential from (1).

The most stringent upper bound on $`\lambda _{231}`$ comes from the leptonic decay width of the tau lepton ,

λ231

<

0.07meR~100GeV,

<

subscript𝜆2310.07subscript𝑚~subscript𝑒𝑅100𝐺𝑒𝑉\lambda_{231}\mathrel{\vbox{\hbox{$<$}\nointerlineskip\hbox{$\sim$}}}0.07{\frac{m_{\widetilde{e_{R}}}}{100GeV}}, (3)
since the selectron exchange diagram through the first and second terms in (2) contributes to the process $`\tau \nu _\mu \nu _\tau \mu `$. The upper limit is not so severe if we compare above bound (3) with the bounds for other RB couplings, e.g., $`\lambda _{133}`$

<

<

\mathrel{\vbox{\hbox{$<$}\nointerlineskip\hbox{$\sim$}}} $`0.0060`$ $`\sqrt{m_{\stackrel{~}{\tau }}/100\mathrm{G}\mathrm{e}\mathrm{V}}`$ derived from the experimental limit on the electron neutrino mass . This is the reason why we consider the $`\widehat{L}_2\widehat{L}_3\widehat{E}_1`$ superpotential.

We are aware of constraints on the slepton masses from the precision measurements at LEP2 . The present lower mass limits at $`95\%`$ C.L. are mν~μ,τ

>

65

>

subscript𝑚subscript~𝜈𝜇𝜏65m_{\widetilde{\nu}_{\mu,\tau}}\mathrel{\vbox{\hbox{$>$}\nointerlineskip\hbox{$\sim$}}}65GeV and me~R

>

69

>

subscript𝑚subscript~𝑒𝑅69m_{\widetilde{e}_{R}}\mathrel{\vbox{\hbox{$>$}\nointerlineskip\hbox{$\sim$}}}69GeV, which are obtained by the analyses on the pair production of the sleptons with the RB decay modes. In the following analysis we take typical input parameters $`m_{\stackrel{~}{\nu }}=100`$GeV, $`\lambda _{231}=0.1`$ and additionally we assume $`Br(\stackrel{~}{\nu }_\tau e\mu )=1`$ for simplicity. This assumption corresponds to the case of a large SU(2) gaugino mass $`M_2`$ and a large SUSY Higgs mass $`\mu `$, e.g., $`M_2>300`$GeV and $`\mu >100`$GeV.

## 3 Production processes

We consider the single sneutrino production as a signal process,

$$e^{}p\mu ^{}\stackrel{~}{\nu }_\tau X.$$
(4)
Note that the mass threshold of the process should be seq

>

mν~

>

subscript𝑠𝑒𝑞subscript𝑚~𝜈\sqrt{s_{eq}}\mathrel{\vbox{\hbox{$>$}\nointerlineskip\hbox{$\sim$}}}m_{\widetilde{\nu}} and it is lower than production threshold at $`e^+e^{}`$ or $`pp`$ colliders as far as the $`\stackrel{~}{\nu }`$ pair production is concerned. The $`\stackrel{~}{\nu }_\tau `$ decays into $`e^{}\mu ^+`$ via the R-parity violating interaction (2). Then we have the final state ***The same final state can be expected when we take the non-zero $`\lambda _{121}`$ or $`\lambda _{122}`$ couplings. For this case $`\stackrel{~}{\nu }_e`$ or $`\stackrel{~}{\nu }_\mu `$ is singly produced with $`\mu ^{}`$ in the $`eq`$ collisions. However, the magnitudes of $`\lambda _{121}`$ and $`\lambda _{122}`$ are more severely constrained from the charged current universality than $`\lambda _{231}`$

$$e^{}p\mu ^{}\mu ^+e^{}X.$$
(5)
Throughout the present work whole calculations of decay widths and cross sections have been performed by using the GRACE/SUSY system, an automatic computation program for SUSY processes . While the GRACE/SUSY system is originally designed to treat such elementary subprocesses as $`e^+e^{}`$, $`eq`$ and $`qq`$ collisions, we have recently succeeded in extending it to $`ep`$ and $`pp`$ collisions. We use the extended new versions which include an interface to the PDFLIB too. For the parton distribution function we have used CTEQ4M . The Feynman diagrams of the subprocesses for(4) as well as relevant SM background are shown in Figs.1 and 2, respectively.

## 4 Numerical results

We obtain the total cross section of 10fb for the signal process (4), where we take $`\lambda _{231}=0.1`$, $`m_{\stackrel{~}{\nu }_\tau }=100`$GeV. It is too small cross section for us to extract from the SM background with the same final states $`\sigma _{SM}(e^{}p\mu ^{}\mu ^+e^{}X)0.3`$nb. Consequently we should suppress the huge background to extract the signal.

In Fig.3 we show the various transverse momentum ($`P_T(e^{})`$, $`P_T(\mu ^{})`$ and $`P_T(\mu ^+)`$) distributions. We find that the signal cross section is much smaller than the SM background in the $`P_T(\mu ^{})`$ distribution in a whole kinematical region. The main contribution to the SM background comes from the two photon processes (graphs 9, 13 in Fig.2)and the virtual photon emission processes (graphs 1, 5, 17, 21 in Fig.2). This is reflected in the rapid increase in the $`P_T`$ distribution with decreasing $`P_T`$ as shown in Fig.3. On the other hand, in the $`P_T(e^{})`$, $`P_T(\mu ^+)`$ distributions, the signal differential cross section can be comparable to the background in a specific kinematical region, $`P_T(e^{})`$, $`P_T(\mu ^+)`$ $``$ $`50`$GeV. This is due to the fact that the $`e^{}`$ and $`\mu ^+`$ are originated from the 2-body decay of the sneutrino with mass $`100`$GeV. We can see the similar small Jacobian peak from the $`Z`$-boson decay for the SM histgrams in the $`P_T(\mu ^+)`$ and $`P_T(\mu ^{})`$ ($`46`$GeV) distributions.

The Monte Carlo event simulation in $`(P_T(\mu ^+),P_T(e^{}))`$ plane is displayed in Fig.4 assuming the integrated luminosity of $`1`$fb-1. While the number of signal events is small compared to the backgrounds, they are characterized by the large transverse momentum. Based on the above observation, we expect that the background can be suppressed if we impose appropriate kinematical cuts on the final leptons. Specifically, we find following condition is suitable for the purpose,

$`P_T(\mu ^+)>20\mathrm{G}\mathrm{e}\mathrm{V},`$ (6)
$`P_T(e^{})>20\mathrm{G}\mathrm{e}\mathrm{V}.`$ (7)
When we apply the $`P_T`$ cuts, we find that the signal can be seen as an event excess in the relevant invariant mass distributions. In Fig.5 we show $`M(\mu ^+\mu ^{})`$ and $`M(e^{}\mu ^+)`$ distributions with the $`P_T`$ cuts (6) and (7). The almost uniform excess for the $`M(\mu ^+\mu ^{})`$ and a sharp peak at $`M(e^{}\mu ^+)=m_{\stackrel{~}{\nu }_\tau }=100`$GeV can be expected.

## 5 Summary and conclusion

We have investigated the single scalar neutrino production at HERA in the framework of the R-parity breaking SUSY model with $`\lambda _{231}0`$. The signal of the process should be multilepton final states $`e^{}\mu ^+\mu ^{}`$ without missing energies.

The main background would be the basic QED reactions, i.e., two photon processes and the virtual photon emission from initial and final fermions. We find that the background can be suppressed by making the appropriate $`P_T`$ cut for for final $`e^{}`$ and $`\mu ^+`$. Then the signal can be clearly seen as a sharp peak at $`M(e^{}\mu ^+)=m_{\stackrel{~}{\nu }_\tau }`$ in the invariant mass distribution.

We conclude that if the R-parity violating coupling constant is λ231

>

0.1

>

subscript𝜆2310.1\lambda_{231}\mathrel{\vbox{\hbox{$>$}\nointerlineskip\hbox{$\sim$}}}0.1 and the mass of scalar neutrino is mν~τ

>

100

>

subscript𝑚subscript~𝜈𝜏100m_{\widetilde{\nu}_{\tau}}\mathrel{\vbox{\hbox{$>$}\nointerlineskip\hbox{$\sim$}}}100GeV, the HERA could be efficient for exploring the $`\widehat{L}\widehat{L}\widehat{E}`$ type interactions in the R-parity breaking SUSY model.

Acknowledgements

The work was supported in part by the Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology of Japan, No. 10440080 and No. 11206203.