text:
# The Clarens Web Services Architecture

## I Introduction

The ascendance of Grid-like grid\_phys\_today  technologies has been all but necessitated by the sheer volume of data produced in both science and commerce. In response to this increased uptake in High Energy Physics in particular, the traditional batch-oriented implementations globus\_toolkit sun\_grid  using home-grown protocols have started to adapt to an industry-wide move to standardized interfaces and protocols globus\_ogsa .

In this context the CAIGEE caigee  project was started to develop a specific application of Grid technologies to the area of interactive analysis by end-user physicists. A diagram showing Clarens as the interface between distributed clients and Grid services is shown in Figure 1.

In keeping with its public standards focus, clients can include anything from an end-user analysis package like ROOT root , web browsers or even web-enabled PDAs. Other Clarens servers may also act as clients.

Back-end services planned or implemented include data movement, Grid-wide execution planning and scheduling, cluster job scheduling as well as metadata catalogs.

## II Infrastructure

In order to save development time and improve scalability, The Clarens server is implemented as an extension to the Apache apache  multi-process web server using the mod\_python extension in the Python byte-code compiled language. Most CMS sites have these components already installed on processing cluster head nodes as part of the Redhat 7.3-based OS used. Clarens itself is both architecture and platform-dependent by virtue of using Python as an implementation language.

The Clarens architecture is depicted in Figure 2 in the order that requests are processed by the server. Firstly (top), the Apache server receives an HTTP POST or GET request from the client, and invokes Clarens based on the form of the URL specified by the client. Other URLs are handled as usual by the server according to its configuration. Secure Sockets Layer (SSL) encrypted connections are handled transparently by the Apache server, with no special coding needed in Clarens itself to decrypt (encrypt) requests (responses). Encryption of network traffic is optional, however, in cases where it is not required, without exposing client or server credentials.

After the request has been processed, a response is sent back to the client, which is usually encoded as an RPC response, but may also be in the form of binary data. GET requests returns a file or an XML-encoded error message to the client, while XML-RPC xmlrpc  or SOAPsoap -encoded POST requests return a similarly encoded response or error message.

## III Authentication

In keeping with the Grid tradition, Clarens uses a so-called Public Key Infrastructure (PKI) based authentication system that relies on X509 formatted certificates issued by a Certification Authority (CA) along with asymmetric encryption using public and private keys. The authentication protocol of the server is implemented at the application level, thus eliminating the need for a custom security layer on the client side.

If an SSL-encrypted connection is used the client’s certificate is provided to the server as part of the connection negotiation stage. The Apache server passes this information to the Clarens layer. This is the default for browser-based clients which in general have a well developed client-side PKI security infrastructure. In this case the authentication step in initiated by calling the RPC method system.auth2() with no arguments. A user session ID is requested by the client by setting the clarens\_username cookie value in the message header to the requested session ID, and setting the clarens\_password cookie value to BROWSER. Clarens responds by returning its own certificate as well as it’s part of the session ID encoded as an RPC response. In subsequent requests the client must set the clarens\_password to this server session ID. Browsers will automatically send these cookie values in subsequent requests.

In the case of an unencrypted connection, or a client not able to send it’s certificate as part of the connection negotiation phase, the session ID and client certificate must be sent using as the username and password in the the HTTP basic authentication header invoking the system.auth() method. The server responds with a list of (1) its certificate, (2) the server session ID encrypted using the user’s public key, and (3) the client session ID encrypted using the server’s private key. This ensures that only someone in possession of the client’s private key can discover the server session ID, and also that the server is in possession of a private key matching the certificate sent as (1) above. In subsequent requests the client should set the client session ID as username, and server session ID as password in the HTTP basic authentication header.

Once this certificate and session ID exchange is completed, both the client and server certificates can be verified against the publicly available CA certificate chain, knowing that the other party is in possession of a matching private key.

## IV Authorization

Authorization of clients to access server resources (invoking methods, accessing files etc.) is done within the framework of a hierarchical Virtual Organization (VO) with members uniquely identified by their Distinguished Names (DNs) issued by the CAs as part of all X509 certificates.

### IV.1 Virtual Organization

Each Clarens server instance manages a tree-like VO structure, as shown in Figure 3, rooted in a list of administrators. This group, named admins, is populated statically from values provided in the server configuration file on each server restart. The list of group members is cached in a database berkeley\_db , as is all VO information. The admins group is authorized to create and delete groups at all levels.

Each group consists of two lists of DNs for the group members and administrators respectively. Group administrators are authorized to add and delete group members, as well as groups at lower levels. The group structure is hierarchical because group members of higher level groups are automatically members of lower level groups in the same branch.

The example in Figure 3 demonstrates the top-level groups A, B, and C, with second level groups A.1, A.2, and A.3.

A more concrete example might be to define groups CMS, Atlas, LHCb, and Alice, then for CMS, to define CMS.USA, CMS.CERN, CMS.UK, CMS.Germany. At the third level, one might define CMS.USA.Caltech, CMS.USA.UFL, CMS.USA.FNAL. Management for the latter three groups may then be delegated to the institutes themselves, thereby implementing a distributed trust model that has lower maintenance overhead as well as being more representative of the real organizational structure.

As a further optimization, the hierarchical information in the DNs may also be used to define membership, so that only the initial significant part of the DN need to be specified. DNs are structured to include information on the country (C), state/province (ST), locality/city (L), organization (O), organizational unit (OU), common name (CN), and e-mail address (Email). An example DN issued by the DOE Science Grid CA is /O=doesciencegrid.org/OU=People/CN=John Smith 12345 for individuals and /O=doesciencegrid.org/OU=Services/CN=host /www.mysite.edu for servers. To add all individuals to a particular group, only /O=doesciencegrid.org/OU=People need to specified as a member DN.

### IV.2 Access Control Lists

Execution of RPCs (web service methods) as well as mapping of certificate DNs to users on the server system is controlled by a set of hierarchical access control lists (ACLs) in a similar fashion to the VO structure described above, and modeled after the access control (.htaccess) files used by Apache.

Methods have a natural hierarchical structure modeled after Python’s module infrastructure. In fact all Clarens modules are also Python modules. Clarens places no arbitrary restrictions on the depth of this hierarchy, but a depth of two or three levels is most common, e.g. module.method or module.submodule.method.

An ACL consist of an evaluation order specification (allow, deny or deny, allow) followed by a list of DNs allowed, groups allowed, DNs denied and groups denied access. A DN or group granted access to a higher level method automatically has access to a lower level method, unless specifically denied at the lower level, and is denied at the higher level unless allowed at a lower level. The ACL specification is therefore evaluated from the lowest applicable level to the highest.

According to the example ACLs in Table 1, when the method mod.meth is invoked, the second ACL is applied, but when any other method in module mod is invoked, the first ACL is applied.

### IV.3 User Mapping

In the traditional batch-oriented Grid architecture, being able to start long-running CPU intensive jobs as a certain system user is of crucial importance. E.g. the Globus toolkit globus\_toolkit  implements the concept of a so-called gridmap file that maps system users to certificate DNs. It is implemented as a flat text file with two values per line, namely the DN and the system user.

Clarens similarly contains the notion of mapping DNs to user names. Instead of a flat file, a structure similar to the ACLs described above is used, with the username taking the place of the method name, with one exception: the deny group and deny individual fields are not used, since denying access to process creation methods can be done using the method ACLs themselves.

At this point it should be clear that in both the VO, ACL, and user mappings specification as little information as possible is stored in order to minimize search times for list memberships. Searching these lists are in the critical path for the invocation of any method, though, and must be optimize as far as practical. This is done by storing lists of DNs as strings in a structure called a ternary tree which is discussed in the appendix.

### IV.4 Auditing

All method invocations are logged with a time-stamp in the Apache server log files to provide a record of client/server transactions.

## V Modular Architecture

Clarens provides a framework for extending its functionality via new modules installed in subdirectories. Each subdirectory would appear as a new method, root, e.g. a system subdirectory would have its methods accessible as system.method etc. Modules can be implemented as either interpreted Python bytecode, or as compiled C/C++ shared libraries if code execution speed is important, and are loaded on demand by the Python interpreter in each Apache process, thereby providing crash-protection between different processes, a major consideration for highly available servers.

Provision is also made for users on the server system to install their own modules in subdirectories under their home directories. These modules are accessible using the format $``$user.module.method.

Figure 4 shows an example of the simplest module implementation for a method echo, that when placed in a file named echo/\__init\__.py returns its argument unchanged when called as echo.echo(argument).

Two important elements identifies this as a Clarens module: the methods\_list and methods\_sig structures identifying the list of methods and their signatures. Work is underway to use WSDL as method signatures instead if the more limited XML-RPC signatures used in this example.

## VI Persistency

Since the HTTP protocol does not require111 In the HTTP 1.1 standard, persistent connections are the default for performance reasons, but the protocol is still inherently stateless, as opposed to e.g. the FTP protocol. persistent connections, it is important that session information be stored persistently on the server side. This has the positive side-effect of allowing clients to survive server failures or restarts transparently without having to re-authenticate themselves to the server in those cases.

The most important session information that is stored by the Clarens server is the authentication information for each session. The Berkeley database berkeley\_db  is used for this purpose.

## VII Scalability, failover, clustering

Since Clarens is built upon commodity software components and standard operating system services, it relies on these components to be set up to achieve these goals. Specifically, the Apache web, server, the filesystem, database and network components needs to be configured by the system administrator.

## VIII Future developments

Work is underway to extend Clarens from being an essentially client/server system to being a truly distributed system in a network of mutually aware peers and superpeers that provide services. The most pressing need for large scale Physics analysis is the need for truly distributed data catalogs for a variety of Physics data and metadata, and of course a matching search capability.

A switch to a relational database for persistent data storage is planned to support more advanced data management than the Berkeley database’s key/value mechanism can provide. Along with this change, file ACLs will also be implemented.

Werever practical Clarens aims to be compatible with the OGSA framework, with support for SOAP and WSDL being important first steps in that direction.

## IX Conclusion

Clarens is a powerful, yet simple web services architecture with a strong emphasis on security for distributed virtual organizations. It draws upon a rich base of commodity protocols and software components to provide a platform for the deployment of analysis-oriented web services.

A companion paper describes services and clients that take advantage of this platform clarens\_apps .

###### Acknowledgements.

This work supported by Department of Energy contract DE-FC02-01ER25459, as part of the Particle Physics DataGrid project ppdg , and under National Science Foundation Grant No. 0218937. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors, and do not necessarily reflect the views of the National Science Foundation. Clarens development is hosted by SourceForge.net sourceforge .

## Appendix: Ternary trees

As pointed out in Section IV, several lists must be searched for every method invocation to enforce access control, establish group membership and map certificate DNs to system user names.

The keys for these lists (method, group, or user names) are kept in Python dictionary objects, which are implemented as hash tables. Although it is certainly possible to store the aforementioned lists in hash tables, these structures lack an important feature, namely the ability to search for initial substrings. E.g. if we are presented with the query DN /O=myorg/OU=People/CN=John Smith, but only the string /O=myorg/OU=People is stored.

Initial substring searches are particularly suited to a data structure called a ternary tree ternary\_tree  that contains nodes for less than, equals and greater than comparisons for string fragments. In contrast to the reference implementation, the Clarens ternary trees contain two characters per node, reducing the storage requirements of the tree.

When presented with the above DN, the tree can be traversed by comparing and branching every two characters, and returning the number of characters matched if a leaf node is reached, or a failure can be signaled if a leaf node is not reached. The returned number can be compared with the length of the query DN: an exact string match corresponds to an equal number of matched characters, while a lower number corresponds to a substring match.

The ternary tree structure is implemented as a Python extension module in C. Informal tests with 10,000 DNs on a 933 MHz PC produces roughly 300,000 searches/second for the worst case where each DN is present in the tree, i.e. the maximum number of branches taken per search.


format: text

text:
# Nonholonomic constraints and Voronec’s equations

## Abstract

Is it allowed, in the context of the Lagrange multiplier formalism, to assume that nonholonomic constraints are already in effect while setting up Lagrange’s function? This procedure is successfully applied in a recent book \[L. N. Hand and J. D. Finch, Analytical Mechanics\] to the problem of the rolling penny, but it does not work in general, as we show by means of a counterexample. It turns out that in many cases the use of nonholonomic constraints in the process of construction of the Lagrangian is allowed, but then the correct equations of motion are the little known Voronec’s equations.

In the Lagrange multiplier formalism as applied to nonholonomic systems, the Lagrangian is written as if there were no constraints. The nonholonomic constraints are taken into account in the formulation of the equations of motion, but not during the construction of the Lagrangian. Setting up the Lagrangian assuming that the constraints are already in effect is completely equivalent to substituting the constraint equations into the Lagrangian written as if there were no constraints. It is tempting to take it for granted that the ensuing reduced Lagrangian together with the relevant constraint equations always lead to the correct equations of motion for the system. The procedure just described is successfully used in a recent book  to solve the problem of a penny rolling on an inclined plane. Unfortunately, contrary to what the mentioned book appears to suggest, this approach is not valid in general, as we proceed to show with the help of a counterexample.

Consider a homogeneous sphere rolling without slipping on a horizontal plane. This problem is treated by the Lagrange multiplier method in . Let $`X,Y,Z`$ be cartesian axes fixed in space with the $`Z`$-axis perpendicular to the plane. The principal moments of inertia with respect to the center of the sphere are all equal to $`2mR^2/5`$. With $`x,y`$ the coordinates of the center of the sphere, the Lagrangian, being equal to the kinetic energy, is given by

$$L=\frac{m}{2}(\dot{x}^2+\dot{y}^2)+\frac{mR^2}{5}𝝎^2.$$
(1)
The constraint equations are

$$\dot{x}=R\omega _y=R\dot{\theta }\mathrm{sin}\varphi R\dot{\psi }\mathrm{sin}\theta \mathrm{cos}\varphi ,\dot{y}=R\omega _x=R\dot{\theta }\mathrm{cos}\varphi R\dot{\psi }\mathrm{sin}\theta \mathrm{sin}\varphi .$$
(2)
In terms of the Euler angles $`\varphi ,\theta ,\psi `$ the Lagrangian (1) takes the form

$$L=\frac{m}{2}(\dot{x}^2+\dot{y}^2)+\frac{mR^2}{5}(\dot{\varphi }^2+\dot{\theta }^2+\dot{\psi }^2+2\dot{\varphi }\dot{\psi }\mathrm{cos}\theta ).$$
(3)
According to the method employed in , which takes into account the rolling constraint in the construction of the Lagrangian, the kinetic energy is written in terms of the rotational degrees of freedom alone by taking the moments of inertia with respect to the contact point of the sphere with the plane. The Lagrangian is now

$$\overline{L}=\frac{1}{2}\frac{7mR^2}{5}(\omega _x^2+\omega _y^2)+\frac{1}{2}\frac{2mR^2}{5}\omega _z^2.$$
(4)
This is exactly the Lagrangian that one obtains by inserting the constraint equations (2) into the Lagrangian (1), which was written as if there were no constraints. In terms of the Euler angles the reduced Lagrangian (4) becomes

$$\overline{L}=\frac{7mR^2}{10}(\dot{\theta }^2+\dot{\psi }^2\mathrm{sin}^2\theta )+\frac{mR^2}{5}(\dot{\varphi }^2+\dot{\psi }^2\mathrm{cos}^2\theta +2\dot{\varphi }\dot{\psi }\mathrm{cos}\theta ).$$
(5)
Since the variables $`x`$ and $`y`$ do not appear in $`\overline{L}`$, according to the reasoning in  the constraint equations (2) are no longer relevant to the formulation of the equations of motion. In particular, the Lagrange equation for $`\theta `$ is

$$\frac{d}{dt}\left(\frac{\overline{L}}{\dot{\theta }}\right)\frac{\overline{L}}{\theta }=0\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}7}\ddot{\theta }5\dot{\psi }^2\mathrm{sin}\theta \mathrm{cos}\theta +2\dot{\varphi }\dot{\psi }\mathrm{sin}\theta =0.$$
(6)
The treatment of this problem by the Lagrange multiplier method shows that the two Lagrange multipliers vanish . According to equation (58$`d`$) of  the correct equation of motion for $`\theta `$ is

$$\ddot{\theta }+\dot{\varphi }\dot{\psi }\mathrm{sin}\theta =0,$$
(7)
which is completely different from equation (6). The differential equations (6) and (7) generally yield different solutions for $`\theta `$ because $`\varphi ,\theta ,\psi ,\dot{\varphi },\dot{\theta },\dot{\psi }`$ can be arbitrarily chosen at any particular instant $`t_0`$. It is clear, therefore, that the approach suggested in  lacks generality, since it works for the rolling penny but fails for the rolling sphere.

In a previous paper  we remarked that it is possible to perform a reduction of the Lagrangian taking into account the constraints, but in this case the correct equations of motion are Voronec’s equations. Given a dynamical system described by the configuration variables $`q_1,\mathrm{},q_n`$, suppose the first $`m`$ velocities are independent and the $`k=nm`$ remaining velocities can be expressed in terms of the independent ones by means of the equations

$$\dot{q}_{m+l}\underset{j=1}{\overset{m}{}}a_{lj}\dot{q}_j=0,l=1,\mathrm{},k,$$
(8)
where the coefficients $`a_{lj}`$ are functions of the generalized coordinates $`q_1,\mathrm{},q_n`$.

Let $`L`$ be the Lagrangian written without taking into account the nonholonomic constraint equations (8). If the last $`k`$ velocities are eliminated from the Lagrangian by means of equations (8), a reduced Lagrangian $`\overline{L}`$ results:

$$L(q_1,\mathrm{},q_n,\dot{q}_1,\mathrm{},\dot{q}_n,t)=\overline{L}(q_1,\mathrm{},q_n,\dot{q}_1,\mathrm{},\dot{q}_m,t).$$
(9)
Voronec’s equations of motion are

$$\frac{d}{dt}\left(\frac{\overline{L}}{\dot{q}_i}\right)\frac{\overline{L}}{q_i}=\underset{\nu =1}{\overset{k}{}}\frac{\overline{L}}{q_{m+\nu }}a_{\nu i}+\underset{\nu =1}{\overset{k}{}}\underset{j=1}{\overset{m}{}}\frac{L}{\dot{q}_{m+\nu }}b_{ij}^\nu \dot{q}_j,i=1,\mathrm{},m,$$
(10)
where

$$b_{ij}^\nu =\frac{a_{\nu i}}{q_j}\frac{a_{\nu j}}{q_i}+\underset{\mu =1}{\overset{k}{}}\left(\frac{a_{\nu i}}{q_{m+\mu }}a_{\mu j}\frac{a_{\nu j}}{q_{m+\mu }}a_{\mu i}\right).$$
(11)
Setting $`q_1=\varphi ,q_2=\theta ,q_3=\psi ,q_4=x,q_5=y`$, in the present case $`m=3`$ and $`k=2`$. The constraint equations (2) can be written in the form (8) with

$$a_{11}=0,a_{12}=R\mathrm{sin}\varphi ,a_{13}=R\mathrm{sin}\theta \mathrm{cos}\varphi ,a_{21}=0,a_{22}=R\mathrm{cos}\varphi ,a_{23}=R\mathrm{sin}\theta \mathrm{sin}\varphi .$$
(12)
The definition (11) furnishes immediately the only nonvanishing coefficients $`b_{ij}^\nu `$:

$$b_{12}^1=b_{21}^1=R\mathrm{cos}\varphi ,b_{13}^1=b_{31}^1=R\mathrm{sin}\theta \mathrm{sin}\varphi ,b_{23}^1=b_{32}^1=R\mathrm{cos}\theta \mathrm{cos}\varphi ,$$
(13)
$$b_{12}^2=b_{21}^2=R\mathrm{sin}\varphi ,b_{13}^2=b_{31}^2=R\mathrm{sin}\theta \mathrm{cos}\varphi ,b_{23}^2=b_{32}^2=R\mathrm{cos}\theta \mathrm{sin}\varphi .$$
(14)
It follows that Voronec’s equation for $`\theta `$ is

$$\frac{7mR^2}{5}\ddot{\theta }mR^2\dot{\psi }^2\mathrm{sin}\theta \mathrm{cos}\theta +\frac{2mR^2}{5}\dot{\varphi }\dot{\psi }\mathrm{sin}\theta =m\dot{x}(b_{21}^1\dot{\varphi }+b_{23}^1\dot{\psi })+m\dot{y}(b_{21}^2\dot{\varphi }+b_{23}^2\dot{\psi }).$$
(15)
The use of the constraint equations (2) and a little algebra reduce the above equation to

$$\ddot{\theta }+\dot{\varphi }\dot{\psi }\mathrm{sin}\theta =0,$$
(16)
which coincides with the correct equation for $`\theta `$ furnished by the Lagrange multiplier method. It is equally straightforward to check that the remaining Voronec equations for $`\varphi `$ and $`\psi `$ coincide with those obtained by the Lagrange multiplier method.

In short, the method employed in  is not valid in general and should not be taught to students. The correct result obtained for the rolling penny is the product of a mere accident.

As a general rule, the use of nonholonomic constraints while setting up the Lagrangian is allowed, but the correct equations of motion are Voronec’s equations, and not the ones given by the Lagrange multiplier method.


text:
# Global Spread of Infectious Diseases

## 1 Introduction and linear model

Recent outbreaks of atypical pneumonia (SARS)  as well as likely future epidemics and pandemics of influenza  provide ample motivation for the study of the global propagation of infectious diseases. (For an overview of modelling of infectious diseases in humans, see ). In this letter we construct a class of simple models of this phenomena, with particular attention to human mobility and the variation of public health capabilities across national and regional boundaries. Spatial dynamics involving the diffusion equation, of the type studied in  (as in, e.g., the spread of rabies in fox populations), is of less interest to us than regional or national variation of parameters. We treat each region as spatially homogenous and focus on mobility (e.g., via air travel) as the mechanism by which disease is transmitted.

The population variables in our model: $`S`$ = susceptibles, $`I`$ = infected, $`R`$ = recovered and $`D`$ = deceased, are functions of time $`t`$ and of discrete geographical location labelled by an index $`i=1,\mathrm{},N`$. Thus our equations describe the time evolution of a $`4N`$ dimensional vector:

$$\stackrel{}{x}=\left(\begin{array}{c}S_1\\ I_1\\ R_1\\ D_1\\ S_2\\ I_2\\ .\\ .\\ .\end{array}\right)$$
(1)
In the simplest linear model, the dynamics of $`\stackrel{}{x}`$ are characterized by a $`4N\times 4N`$ matrix whose entries are probabilities per unit time: $`t_i`$ = probability of transmission of disease from an infected to susceptible individual in region $`i`$, $`r_i`$ = probability of recovery of infected individual in region $`i`$, $`d_i`$ = probability of death of infected individual in region $`i`$, $`m_{ij}`$ = probability of movement of an infected individual from country $`i`$ to country $`j`$. Each of these probabilities incorporates several additional parameters that would appear in a more complicated model. For example, $`t_i`$ depends on other factors such as the probability and effectiveness of quarantine, the length of time an infected individual remains asymptomatic, the population density, etc. A very simple modification would be to include asymptomatic as well as symptomatic carriers of the infection by expanding our matrix to $`5N\times 5N`$ dimensions. Contagious asymptomatics would have a large transmission probability, as they are difficult to quarantine. We will return later to enhancements of the basic model.

The basic parameters $`r_i,t_i,d_i,m_{ij}`$ vary from region to region, and can be estimated from epidemiological data. Our organizing principle for dividing the world into regions is the relative homogeneity of these parameters within each region. This division may or may not follow national boundaries, as the parameters may vary significantly (e.g., from urban to rural regions) within a given country.

The basic reproductive rate , commonly denoted $`_0`$, is defined as the average number of secondary cases caused by a single infected individual. A simple computation (neglecting migration) for region $`i`$ yields

$$_0^i=t_i+(1r_id_i)t_i+(1r_id_i)^2t_i+\mathrm{}=\left(\frac{t_i}{r_i+d_i}\right).$$
(2)
When $`_0>1`$ the infected population grows exponentially.

The system of equations governing the time evolution of $`\stackrel{}{x}`$ is

$$\frac{d\stackrel{}{x}}{dt}=M\stackrel{}{x}+\stackrel{}{J},$$
(3)
where $`\stackrel{}{J}`$ is a possible source of infected individuals coming from zoonosis, an animal reservoir of disease. All elements of $`\stackrel{}{J}`$ are zero except the $`(4i2)`$ and $`(4i3)`$ entries corresponding to an animal reservoir in region $`i`$ which contributes to $`dI_i/dt`$ and $`dS_i/dt`$. Note that, unlike the well-known SIR model , equation (3) is linear. The transmission rate per unit time is given by $`t_iI_i`$ rather than $`t_iI_iS_i`$. This makes our model slightly less realistic, as the transmission rate per infected individual does not drop as the number of susceptibles goes to zero. The SIR model predicts some fraction of uninfected susceptibles $`S(t\mathrm{})`$ even when $`_0>1`$, whereas our model would predict that each individual in the entire population is eventually be infected. This leads to an overestimate relative to the SIR model of the total number of infected in countries with $`_0>1`$. However, the difference is only numerically significant when $`_0`$ is close to 1. (We do not expect the input parameters to be sufficiently well-determined that, e.g., a factor of two in the predicted number of infections is very significant.) Another method for cutting off the exponential growth is to add (nonlinear) logistic terms of the form $`k_iI_i^2`$ to the right hand side of the equation, at the cost of adding additional parameters.

Equation (3) can be solved analytically:

$$\stackrel{}{x}=M^1J+e^{Mt}\left(\stackrel{}{x}_0+M^1J\right),$$
(4)
where $`\stackrel{}{x}_0`$ is the initial condition. The importance of the variables $`S_i,R_i,D_i`$ is to impose the overall conservation of humans, which is a consequence of the form of $`M`$, and the constraint that $`S,I,R,D`$ be positive semi-definite at all times. Although equations (3) and (4) appear simple, the constraints lead to a nontrivial system.

Below we exhibit the matrix M for the case of two countries, $`i=1,2`$.

$$M=\left(\begin{array}{cccccccc}0& t_1& 0& 0& 0& 0& 0& 0\\ 0& r_1d_1m_{12}+t_1& 0& 0& 0& m_{21}& 0& 0\\ 0& r_1& 0& 0& 0& 0& 0& 0\\ 0& d_1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& t_2& 0& 0\\ 0& m_{12}& 0& 0& 0& r_2d_2m_{21}+t_2& 0& 0\\ 0& 0& 0& 0& 0& r_2& 0& 0\\ 0& 0& 0& 0& 0& d_2& 0& 0\end{array}\right).$$
(5)
Note the sum rule implicit in this matrix: $`_i^{4N}x_i=\mathrm{constant}`$, which reflects conservation of humans.

Despite the simplicity of these models, they may be useful tools for public health policy. For example, the economic impacts on developed country economies from diseases in developing nations can be estimated. This data can be used in cost-benefit assessments of improvements in public health infrastructure in both developed and developing nations.

## 2 Linear model: analytics and simulations

We now discuss some of the analytic properties of our linear models. First, we note that the system of $`4N`$ equations can be reduced to a smaller number of equations, due to the large number of zeroes in $`M`$. In essence, the dynamics involves only the infected individuals $`I_i`$: the evolution of the $`S_i,R_i`$ and $`D_i`$ variables are dependent on $`I_i`$. Thus there are really only $`N`$ equations (but subject to constraints; see below), as can be seen from the fact that $`M`$ has only $`N`$ non-zero eigenvalues. The matrix $`M`$, restricted to the $`N\times N`$ subspace of infected populations $`I_i`$, is diagonal up to the (small) off-diagonal mobility entries $`m_{ij}`$. If the mobility matrix were symmetric, we could diagonalize $`M`$ via an orthogonal change of basis. However, in general $`m_{ij}m_{ji}`$, so we need to use a bi-linear transformation. That is, we can obtain a diagonal matrix $`M_D`$

$$M_D=LMR,$$
(6)
where $`L=R^1=R^T`$ up to corrections which vanish for symmetric $`m_{ij}`$. In the basis $`\stackrel{}{y}=R^1\stackrel{}{x}`$, equation (3) becomes

$$LR\frac{d\stackrel{}{y}}{dt}=M_D\stackrel{}{y}+L\stackrel{}{J}.$$
(7)
While the eigenvalues $`\lambda _i`$ tell us a great deal - essentially, what would happen in each region if it were isolated - there is additional information in the $`L,R`$ matrices which relate the original regional basis to the diagonal basis. Below we investigate explicit scenarios involving a developing country $`i`$ with positive eigenvalue $`\lambda _i`$, and a developed country $`j`$ with eigenvalue $`\lambda _j`$ which is either negative, or, if positive, smaller than $`\lambda _i`$. In these scenarios the $`m_{ij}`$ mobility is much more important than the reverse migration $`m_{ji}`$, since the spread of disease is mainly unidirectional. Hence, we could simply assume $`m_{ij}=m_{ji}`$ without changing the results appreciably. In that case, $`L=R^1`$, so equation (7) becomes uncoupled.

In the $`N=2`$ case, the non-zero eigenvalues are

$`\lambda _1`$ $`=`$ $`r_1d_1m_{12}+t_1`$
$`\lambda _2`$ $`=`$ $`r_2d_2m_{21}+t_2.`$ (8)
Note that a positive eigenvalue $`\lambda _i`$ corresponds to $`_0^i>1`$.

Two generic implications can be stated, depending on the sign of $`\lambda _i`$. Here we assume that the entries of the off-diagonal mobility matrix $`m_{ij}`$ are small compared to the other probabilities, so that we can discuss the eigenvalue associated an individual region $`i`$. This is likely to be the case in any realistic scenario where the mobility matrix reflects air travel.

A country with positive eigenvalue $`\lambda _i`$ will experience exponential growth in the number of infecteds $`I_i`$. This growth is only limited by the constraint of total population: eventually roughly the entire susceptible population has been infected and after a long period of time there are only recovered and dead individuals. (This situation is familiar in the case of the common cold or influenza.) There are two important timescales involved: $`\lambda _i^1`$, which determines the rate of exponential growth, and $`|\lambda _it_i|^1`$, ($`\lambda _it_i`$ is negative) which gives the rate of decline of $`I_i`$ after saturation occurs and there are no more susceptibles to be infected ($`S_i=0`$, so effectively the transmission probability $`t_i`$ can be set to zero in the corresponding eigenvalue). The outbreak in a highly susceptible country cannot last much longer than the sum of these timescales.

A country with a negative eigenvalue is capable of suppressing any outbreaks of the disease. Nevertheless, there may be a steady state number of infections

$$I_j\frac{J_j}{|\lambda _j|}+\underset{ij}{}\frac{m_{ij}I_i}{|\lambda _j\lambda _i|}$$
(9)
due to the animal reservoir $`J`$ and/or the migration of infected individuals from other regions. This equation is valid for $`m_{ij}`$ small compared to $`|\lambda _j\lambda _i|`$. (The difference in eigenvalues results from the diagonalization of $`M`$, and is analogous to the result in quantum mechanics for the change in energy eigenstate under a small perturbation to the Hamiltonian.) The most dangerous sources of disease migration are countries $`i`$ with positive eigenvalue $`\lambda _i`$. The total flux of migration is roughly determined by the mobility $`m_{ij}`$, the population $`S_i(0)`$ and the timescales $`\lambda _i^1`$ and $`|\lambda _it_i|^1`$.

Consider a plausible worst case scenario involving a developing country $`i`$ with positive $`\lambda _i`$ and much larger population than developed country $`j`$ which has negative $`\lambda _j`$. The developing country is the main source of infection for the developed country. In this case, we obtain the simple relation

$$I_j\frac{m_{ij}I_i}{|\lambda _j\lambda _i|}$$
(10)
relating the number of infected in the two countries at any particular time. The total number of infections in country $`j`$ (during the entire epidemic) is related to $`I_j(t)`$ in a complicated way, depending on recovery and death rates. However, by solving equation (4) we can obtain a simple expression for the ratio of total number of infections $`N_i`$ in the two countries

$$\frac{N_j}{N_i}=\frac{m_{ij}}{|\lambda _j|}\left(\frac{d_j+r_j}{d_i+r_i}\right),$$
(11)
to leading order in $`m_{ij}`$. The first ratio in this result can be interpreted as the fraction of infected individuals who travel from country $`i`$ to $`j`$ in the timescale $`|\lambda _j|^1`$, which is roughly the “halving” time in country $`j`$. As long as the death and recovery rates in the two countries are not too different, this can be used as a rule of thumb to estimate the ratio of total number of infections. Obviously it behooves country $`j`$ to make $`|\lambda _j|^1`$ as small as possible.

We can also obtain simple expressions for the peak number of infected individuals in both countries (neglecting animal reservoirs):

$$I_i^{\mathrm{peak}}=\frac{\lambda _i}{t_i}\left(P_iI_i(0)\left(1+\frac{t_i}{\lambda _i}\right)\right)\frac{\lambda _i}{t_i}P_i$$
(12)
and

$`I_j^{\mathrm{peak}}`$ $`=`$ $`{\displaystyle \frac{m_{ij}}{|\lambda _i\lambda _j|}}\left({\displaystyle \frac{\lambda _i}{t_i}}\left(P_iI_i(0)\left(1+{\displaystyle \frac{t_i}{\lambda _i}}\right)\right)+I_i(0)\left({\displaystyle \frac{\lambda _i}{I_i(0)t_i}}\left(P_iI_i(0)\left(1+{\displaystyle \frac{t_i}{\lambda _i}}\right)\right)\right)^{\frac{\lambda _j}{\lambda _i}}\right)`$ (13)
$``$ $`{\displaystyle \frac{m_{ij}}{|\lambda _i\lambda _j|}}{\displaystyle \frac{\lambda _i}{t_i}}P_i`$
where $`I_i(0)`$ is the initial number of infecteds, and $`P_i`$ is the total population, in country $`i`$. This gives the maximum capacity required of the medical system during the epidemic. Again, our result is likely to be an overestimate relative to a more realistic SIR or logistic model.

Below we describe the results of two illustrative simulations of the linear model. The parameters are chosen to be realistic average probabilities per week of transmission, recovery, death and international travel. We integrate equation (3) forward numerically, discretizing time in units of weeks, and starting with a random number of infecteds between 1 and 10. Country 1 has population $`10^9`$ and country 2 has population $`10^8`$. Note that the figures display the infected population $`I_i(t)`$ as a function of time, not the total number of individuals who have ever been infected. The latter quantity, although not displayed, agrees precisely with the result given in equation (11). The peak infected results in equations (12) and (13) are also confirmed.

1) One positive and one negative eigenvalue. The epidemic infects the entire population of country 1 (figure (1)), but infects only a small fraction of the population of country 2 (figure (2)). The parameter values used are $`t_1=1;r_1=0.7;d_1=0.2;t_2=0.85;r_2=0.9;d_2=0.05;m_{12}=0.00001;m_{21}=0`$, eigenvalues $`\lambda _1=0.09999`$, $`\lambda _2=0.1`$. The overall results are consistent with our result (11).

2) Two positive eigenvalues. The epidemic sweeping through country 1 (figure (3)) drives rapid growth in country 2 (figures (4),(5)) until it has completed its course. At this point slower growth characterized by $`\lambda _2`$ (which is very small) resumes. In the end all susceptible individuals in both countries are infected. The falloff of $`I_2`$ appears abrupt in figure (5) because $`|\lambda _2t_2|`$ is much larger than $`\lambda _2`$. $`t_1=1;r_1=0.75;d_1=0.2;t_2=0.851;r_2=0.8;d_2=0.05;m_{12}=0.00001;m_{21}=0`$, eigenvalues $`\lambda _1=0.04999`$, $`\lambda _2=0.001`$.

The numerical values of parameters have been varied only slightly between the two examples, but the resulting behaviors are very different, largely because the eigenvalue $`\lambda _2`$ has changed sign.

## 3 Nonlinear models

To make the models more realistic, we can let the parameters $`t_i,r_i,d_i`$ depend on local quantities $`x_i`$. The resulting equations are nonlinear, although still straightforward to simulate numerically. For example, it is possible that a region’s medical capabilities are overwhelmed as the number of infected individuals increases. Indeed, in both the SARS outbreak and the influenza of 1918 (Spanish flu) medical personnel were disproportionately affected.

As a simple example, we take the parameters from simulation 1 in the previous section, and let the recovery rate $`r_2`$ of country 2 (which had a negative eigenvalue) interpolate between the initial value $`r_2`$ when $`I_2=0`$ and $`r_2^{}`$ when $`I_2`$ is large, according to the formula

$$r_2(I_2)=r_2+(r_2^{}r_2)\left(\frac{I_2}{I_2+I_2^{}}\right).$$
(14)
This reflects a saturation of medical treatment capabilities after some characteristic number of infections. When $`I_2`$ is much larger than $`I_2^{}`$ the recovery probability per unit time is reduced from its initial value of $`r_2`$ to $`r_2^{}`$, and the eigenvalue $`\lambda _2`$ may become positive. We take $`I_2^{}=2000`$, $`r_2=.9`$ and $`r_2^{}=.7`$. The results are shown in figures (6) and (7). Figure (7) reveals the onset of the nonlinear behavior when $`I_2`$ is a few thousand (compare to figure (2)).

## 4 Internet viruses and worms

These models can be modified to describe the propagation of computer viruses and worms. (From the point of view of our models, there is no essential difference between the two.) A virus commonly inserts itself into other program files, in the same manner that a virus in nature takes over the workings of normal cells. When the infected program runs, the virus code gets a chance to inspect its environment and look for and infect new carriers in the form of other program files. A worm is a self-replicating program that does not alter files but resides in active memory and duplicates itself by means of computer networks. Worms use facilities of an operating system that are meant to be automatic and invisible to the user. It is common for worms to be noticed only when their uncontrolled replication consumes system resources, slowing or halting other tasks.

The variation of parameters by region in our model would be a consequence of different levels of security administration in home, small business, educational and corporate networks.

Mobility $`m_{ij}`$ has no analog here, as the infected computer does not hop from network to network. Rather, an infected machine can infect other machines connected to the Internet, leading to a non-local (off-diagonal) transmission matrix $`t_{ij}`$. Machines on local networks $`k`$ with strong firewalls and virus scanning would be difficult to infect, leading to correspondingly small entries $`t_{ik}`$. However, once a virus or worm has penetrated the firewall, spread within a local network can be quite rapid, so the diagonal elements $`t_{kk}`$ may be quite large. Often, the only factor limiting the rate of infection is saturation of available bandwidth by worm scanning activity.

“Slammer”, the most successful worm to date, exploited a security hole in Microsoft’s SQL Server (a database server), and infected 75,000 machines around the world in less than a half hour on January 25, 2003 . The initial doubling time of the infected population was 8.5 seconds, and the growth curve displayed a typical logistic shape.

## 5 Discussion

We have attempted to model the spread of infectious diseases in a global environment, taking into account geographical variation of public health capabilities as well as human mobility (i.e., via air travel). We believe our models capture a large range of possible behaviors using a minimal set of parameters, each of which can be roughly estimated from data.

One interesting conclusion from our models is that typical international mobility – the probability per unit time of international travel for a given infected individual, estimated at $`m_{ij}10^5`$ per week – is still sufficiently small that a country with well-developed public health infrastructure (effectively, a negative eigenvalue $`\lambda `$) can resist an epidemic even when other more populous countries experience complete saturation. In the quasi-realistic simulation 1 (figures (1),(2)), of order $`10^5`$ infections occur in country 2, even though the disease has swept completely through country 1. In reaching this conclusion, we kept the mobility parameter fixed during the outbreak, and did not assume any draconian quarantine on international travellers arriving in country 2. Such measures would reduce the number of infections in country 2 considerably. Of course, this conclusion assumes that the public health infrastructure in country 2 remains robust during the outbreak. In the nonlinear simulation 3 (figures (6), (7)), we see that a breakdown in the medical system can lead to grave consequences.

In the case of two countries, one of which is a “reservoir” with positive eigenvalue $`\lambda _i`$ and the other with negative eigenvalue $`\lambda _j`$, a good rule of thumb arising from equation (11) is that the ratio of total number of infections in the two countries is given by the fraction of infected individuals who migrate in a timescale $`|\lambda _j|^1`$, which is the “halving” time for the epidemic in country $`j`$. In our simulation 1, this timescale is about two months, and the fractional mobility over that period is $`10^4`$, leading to $`10^5`$ infections in country 2 if the entire $`10^9`$ population of country 1 is infected. The maximum number of infections at any given time in the simulation is $`510^3`$ (figure (2)). If the medical system of country 2 can treat this number of patients without breaking down (entering the nonlinear regime), it can prevent a larger outbreak.

Our formula (11) is not applicable to early SARS data since the eigenvalues have clearly been time-dependent during the early stages of the epidemic (estimates of $`_0`$ vary widely during different time periods ). Even developed countries like Canada and Taiwan exhibited $`_0>1`$ during the first months (positive eigenvalues), although as of this writing eigenvalues appear to be negative for all countries.

## Acknowledgements

AZ thanks A. Madhav for conversations. SH is supported under Department of Energy contract DE-FG06-85ER40224. AZ is grateful for the warm hospitality of the University of Oregon, where this work was begun, and acknowledges support from NSF grant PHY 99-07949.


text:
# Pendulum View of The Schrödinger Equation

## Appendix A Winding number $`W(\lambda )`$

For a given force function $`A(x)`$, $`W(\lambda ^{})W(\lambda )`$ if $`\lambda ^{}>\lambda `$.

Proof: Assume the opposite, $`W(\lambda ^{})<W(\lambda )`$. Since the difference is at least one, the pendulum must make more rotations at $`\lambda `$ than at $`\lambda ^{}`$. Note also that at $`x\mathrm{}`$ one has $`\alpha (x,\lambda ^{})>\alpha (x,\lambda )`$. These two facts infer that there must exist a point $`\stackrel{~}{x}`$, where the two pendulum motion functions, $`\alpha (x,\lambda ^{})`$ and $`\alpha (x,\lambda )`$ cross each other for the first time. At this crossing position, $`x=\stackrel{~}{x}`$, one should have $`\frac{d\alpha }{dx}(\stackrel{~}{x},\lambda ^{})<\frac{d\alpha }{dx}(\stackrel{~}{x},\lambda )`$. However, this contradicts the fact,

$$\frac{d\alpha }{dx}(\stackrel{~}{x},\lambda ^{})\frac{d\alpha }{dx}(\stackrel{~}{x},\lambda )=2(\lambda ^{}\lambda )>0.$$
(18)
So, the assumption is not true. We have $`W(\lambda ^{})W(\lambda )`$.

## Appendix B Transformations of Zakharov-Shabat equations

The Zakharov-Shabat (ZS) equations are

$`idU_1(x)/dx+\varphi (x)U_2(x)`$ $`=`$ $`\lambda U_1(x),`$ (19)
$`idU_2(x)/dx\varphi ^{}(x)U_1(x)`$ $`=`$ $`\lambda U_2(x),`$ (20)
where the complex function $`\varphi (x)`$ acts as a potential. We choose to write $`\varphi (x)=A(x)e^{iS(x)}`$, where $`A(x)>0`$ and $`S(x)`$ are real.

Schrödinger equation.– Introduce a transformation

$$\psi _\pm (x;\lambda )=\frac{U_1(x)e^{iS(x)/2}\pm iU_2(x)e^{iS(x)/2}}{\sqrt{2\lambda +S^{}(x)}},$$
(21)
where $`S^{}=dS/dx`$. By straightforward computation, one can verify that $`\psi _\pm `$ satisfy the Schrödinger equation

$$d^2\psi _\pm /dx^2+V_\pm (x)\psi _\pm =(\lambda +S^{}/2)^2\psi _\pm ,$$
(22)
where

$$V_\pm (x)=(A^2\pm A^{})\pm \frac{AS^{\prime \prime }}{2\lambda +S^{}}\frac{2S^{(3)}(\lambda +S^{})3(S^{\prime \prime })^2}{4(2\lambda +S^{})^2}.$$
(23)
In the special case $`S(x)=\mathrm{constant}`$, the above equation becomes

$$d^2\psi _\pm /dx^2+V_\pm (x)\psi _\pm =\lambda ^2\psi _\pm ,$$
(24)
with much simpler potentials

$$V_\pm (x)=A^2\pm A^{}.$$
(25)
These two potentials $`V_\pm `$ are $`\stackrel{~}{V}`$ and $`V`$ used in the main text.

Pendulum equation.– For the discrete eigenvalues, the ZS equations have $`|U_1|=|U_2|=\sqrt{\rho }`$. This allows one to introduce another transformation with

$$U_1=\sqrt{\rho }e^{i\theta _1},U_2=\sqrt{\rho }e^{i\theta _2}.$$
(26)
Substituting it into the ZS equations, one obtains

$`d\rho /dx`$ $`=`$ $`2A\rho \mathrm{sin}(\theta S),`$ (27)
$`d\theta /dx`$ $`=`$ $`2\lambda +2A\mathrm{cos}(\theta S),`$ (28)
where $`\theta =\theta _1\theta _2`$. The above two equations become Eqs.(3,4) when $`S(x)=0`$ with $`\alpha =\pi /2\theta `$.

The combination of Eq.(21) and Eq.(26) results in the transformations Eqs.(2,16).


text:
# Galaxy Types in the Sloan Digital Sky Survey Using Supervised Artificial Neural Networks

## 1 Introduction

The comparison of the observed distribution of galaxies and their properties with that predicted by theory is an important task in cosmology. In recent years datasets have become available which enable the comparison to include large samples and detailed galaxy parameters. The Sloan Digital Sky Survey (SDSS, York et al., 2000) provides a dataset of unprecedented size and quality and thus enables significant improvement in the detail of the comparison.

One can measure an almost limitless number of parameters to describe a galaxy. It is desirable to have as much information as possible in the fewest parameters, either continuous or discrete. A one parameter galaxy ‘type’ is particularly convenient. Examples are the well-known Hubble system, or spectral types based on lines or principal component analysis.

Principal component analysis (PCA), Fisher Matrix and other techniques provide a linear method of reducing the dimensionality of the parameter space in this way. However galaxy parameters are in general correlated in non-linear ways, thus a non-linear approach may be more appropriate. Various methods exist, including non-linear PCA (e.g. http://www.cis.hut.fi/projects/ica/), Information Bottleneck (Slonim et al., 2001), and artificial neural networks (ANNs). The latter approach is adopted here.

The derived parameters should be physically meaningful, i.e. they should be directly predicted by theories of galaxy and large scale structure formation, or be related in a quantitative way. For PCA numerous studies have found that the principal components of galaxy spectra correlate with various physical processes such as star formation (via absorption and emission line strengths of, for example, the H$`\alpha `$ line), and to galaxy colour and morphology. PCA has been applied to the SDSS and yields a one parameter spectral type known as the eClass (Connolly & Szalay, 1999). A similar parameterization, the $`\eta `$ class, has been made for the 2dF galaxy redshift survey (Madgwick et al., 2002).

Here ANNs in the Matlab Neural Network Toolbox environment (http://www.mathworks.com/) are used to map galaxy parameters from Data Release One (DR1) of the SDSS on to a single continuous ‘type’. Here we consider three different types: morphological classification, spectral type and redshift, with standard photometric parameters as input.

Previous studies involving galaxy classification using ANNs include Storrie-Lombardi et al. (1992), Serra-Ricart et al. (1993), Adams & Woolley (1994), Lahav et al. (1995), Naim et al. (1995), Folkes, Lahav & Maddox (1996), Lahav et al. (1996), Odewahn et al. (1996), Naim, Ratnatunga & Griffiths (1997a, b), Molinari & Smareglia (1998), de Theije & Katgert (1999), Windhorst et al. (1999), Bazell (2000), Bazell & Aha (2001), Ball (2001), Goderya & Lolling (2002), Odewahn et al. (2002), Cohen et al. (2003) and Madgwick (2003). However none of these used a dataset of the size and quality of DR1, or the Levenberg-Marquardt training algorithm (§3), widely used in neural network research.

The layout of the rest of this paper is as follows: in §2 the SDSS is summarized and the datasets used are described. In §3 we describe the ANNs; §4 presents the results, followed by discussion in §5 and conclusions in §6.

## 2 Data

The SDSS is a project to map $`\pi `$ steradians of the northern galactic cap in five bands ($`u`$, $`g`$, $`r`$, $`i`$ and $`z`$) from 3500–8900 Å. This will provide photometry for of order $`5\times 10^7`$ galaxies (Fukugita et al. 1996; Gunn et al. 1998; Lupton et al. 2001; Hogg et al. 2001; Smith et al. 2002; Pier et al. 2003). A multifibre spectrograph will provide redshifts and spectra for approximately $`10^6`$ of these. A technical summary of the survey is given in York et al. (2000).

The data released to the community so far consists of the June 2001 Early Data Release (EDR, Stoughton et al., 2002) and the April 2003 Data Release 1 (DR1, Abazajian et al., 2003). These respectively provide photometric parameters and images for one and several million galaxies and spectra for 39,959 and 134,015 galaxies. This paper uses galaxies from DR1.

The SDSS galaxies with spectra consist of a ‘main’, flux-limited sample ($`r<17.77`$), with a median redshift of 0.104 (Strauss et al., 2002) and a luminous red galaxy sample, approximately volume-limited to $`z0.4`$ (Eisenstein et al., 2001). Only the main sample galaxies are used here.

### 2.1 Galaxy Samples

We used the main galaxy sample from DR1, with sample cuts of reddening corrected $`r`$-band magnitude $`r<17.77`$, confidence in spectroscopic redshift zConf $`>0.85`$ and spectroscopic object class specClass $`=`$ GALAXY or emission line galaxy GAL\_EM. This gave 104,619 galaxies. For each of the training, test and simulation samples (see §3) galaxies with severely outlying parameters ($`>10\sigma `$ from the mean value for the parameter, generally indicative of a measurement error) were iteratively removed for each parameter in turn. 2,240 were removed in this way, leaving 102,379. See §2.2 for a description of parameters used. Galaxies with outlying target types were similarly removed. The order in which the parameters are presented may affect the number of galaxies removed, but the difference is negligible given the small number of objects affected (almost the same outliers are removed whatever the order of parameter presention). The parameters were individually normalized to zero mean and unit variance for input into the neural network (see below). For eClass and redshift the training samples were evened out by binning the galaxies by target type and removing random galaxies from the most populated bins until the maximum number of galaxies in a bin was twice the mean number. Bins with less than this were unaffected. This culling ensures that the training of the network is not dominated by only a small region of parameter space where there are large numbers of galaxies, which worsens the performance on the rest of the space, and left a total of 98,402 galaxies. The culling does remove the Bayesian prior of the relative number of each type of galaxy, but the training samples are large enough that the performance on the test sample is improved rather than hindered. A more sophisticated method of creating an even sample is to use K-means clustering or a self-organizing map (e.g. Tagliaferri et al., 2002).

### 2.2 Galaxy Parameters

The parameters used as input to the neural networks, all available in DR1, are shown in Table 1.

The magnitudes are corrected for galactic reddening using the corrections derived from Schlegel, Finkbeiner & Davis (1998).

The galaxy images are fitted with the de Vaucouleurs profile (de Vaucouleurs, 1948)

$$I(r)=I_0\mathrm{exp}\{7.67[(r/r_\mathrm{e})^{1/4}1]\},$$
(1)
and the exponential profile (Freeman, 1970)

$$I(r)=I_0\mathrm{exp}(1.68r/r_\mathrm{e}),$$
(2)
where $`I_0`$ and $`I(r)`$ are the intensities at radii 0 and $`r`$, and $`r_\mathrm{e}`$ is the half-light radius for the galaxy. The profiles are truncated to go smoothly to zero at $`8r_\mathrm{e}`$ and $`4r_\mathrm{e}`$ respectively.

The profile likelihoods are standard $`\chi ^2`$ fits. The model magnitude is that from the better of the two fits.

The Petrosian magnitude is a modified form of that introduced by Petrosian (1976). It measures a constant fraction of the total light. The Petrosian flux is given by

$$F_\mathrm{P}_0^{N_\mathrm{P}r_\mathrm{P}}2\pi r^{}𝑑r^{}I(r^{})$$
(3)
where $`r_\mathrm{P}`$ is the Petrosian radius, which is the value at which the Petrosian ratio of surface brightnesses

$$R_\mathrm{P}(r)\frac{_{0.8r}^{1.25r}2\pi r^{}𝑑r^{}I(r^{})/[\pi (1.25^20.8^2)]}{_0^r2\pi r^{}𝑑r^{}I(r^{})/(\pi r^2)}$$
(4)
has a certain value, chosen in the SDSS to be 0.2. The number $`N_\mathrm{P}`$ of Petrosian radii within which the flux is measured is equal to 2 in the SDSS.

The magnitude $`m`$, as with the model magnitude, is then given in asinh units, which are virtually identical to the usual astronomical magnitudes (Pogson, 1856) at high signal to noise but work at low signal to noise and negative flux:

$$m=\frac{2.5}{\mathrm{ln}10}\left[\mathrm{asinh}\left(\frac{f/f_0}{2b}\right)+\mathrm{ln}b\right],$$
(5)
where $`b`$ is a softening parameter. Further details are given in Lupton, Gunn & Szalay (1999) and Stoughton et al. (2002).

The concentration index is $`R_{50}/R_{90}`$ where $`R_{50}`$ and $`R_{90}`$ are the radii within which 50 and 90 per cent of the Petrosian flux is received.

The surface brightness used here is given by

$$\mu =r+5\mathrm{log}(\pi r_\mathrm{P}^2),$$
(6)
$`r_\mathrm{P}`$ being the Petrosian radius in the $`r`$ band.

Parameters other than magnitudes and colours are measured in the $`r`$ band, since this band is used to define the aperture through which Petrosian flux is measured for all five bands. Further details of all the parameters are given on the DR1 webpage (http://www.sdss.org/dr1).

### 2.3 Target Types

The networks were separately trained on the following three targets.

#### 2.3.1 Eyeball Morphological Type

1875 SDSS galaxies have been classified into morphological types by Nakamura et al. (2003). The system used was a modified version of the T-type system (de Vaucouleurs, 1959), with the types being assigned in steps of 0.5 from 0 (early type) to 6 (late type). Unassigned types ($`1`$) and galaxies flagged as being likely to have bad photometry were removed.

The Nakamura et al. catalogue is based on a pre-DR1 version of SDSS data and so their catalogue was matched to DR1 by equatorial coordinates with a tolerance of 0.36 arcsec, so that the number of duplicate matches is negligible. This gave 1399 matches.

#### 2.3.2 eClass

The eClass is a continuous one parameter type assigned from the projection of the first three principal components (PCs) of the ensemble of SDSS galaxy spectra. The locus of points forms an approximately one dimensional curve in the volume of PC1, PC2 and PC3. This is a generalization of the mixing angle $`\varphi `$ in PC1 and PC2

$$\varphi =\mathrm{tan}^1\left(\frac{a_2}{a_1}\right),$$
(7)
where $`a_1`$ and $`a_2`$ are the eigencoefficients of PC1 and PC2.

The range is from approximately $`1`$ (corresponding to late type galaxies) to 0.5 (early type).

The eClass is also robust to missing data in the spectra used for its derivation, and is almost independent of redshift. Further details can be found in Connolly et al. (1995), Connolly & Szalay (1999), and Yip et al. (2002).

#### 2.3.3 Redshift

The redshift is calculated automatically by the SDSS spectroscopic software pipelines (Stoughton et al., 2002), Frieman et al. (in preparation), and has a success rate of almost 100 per cent.

## 3 Artificial Neural Networks

ANNs, as collections of interconnected neurons each able to carry out simple processing were originally conceived as being models of the brain. This is still true, however the networks used here are vastly smaller and simpler and are best described in terms of non-linear extensions of conventional statistical methods.

The supervised ANN takes parameters as input and maps them on to one or more outputs. A set of vectors of parameters, each vector representing a galaxy and corresponding to a desired output, or target, is presented. The network is trained and is then able to assign an output to an unseen parameter vector.

This is achieved by using a training algorithm to minimize a cost function which represents the difference between the actual and desired output. The cost function $`c`$ is commonly of the form

$$c=\frac{1}{N}\underset{k=1}{\overset{N}{}}(o_kt_k)^2,$$
(8)
where $`o_k`$ and $`t_k`$ are the output and target respectively for the $`k`$th of $`N`$ objects.

In general the neurons could be connected in any topology, but a commonly used form is to have an $`a:b_1:b_2:\mathrm{}:b_n:c`$ arrangement, where $`a`$ is the number of input parameters, $`b_{1\mathrm{}n}`$ are the number of neurons in each of $`n`$ one dimensional ‘hidden’ layers and $`c`$ is the number of neurons in the final layer, equal to the number of outputs. Here we have one output, $`c=1`$. Multiple outputs can give Bayesian a posteriori probabilities that the output is of that class given the values of the input parameters. (This is classification, whereas a single output, $`c=1`$, is strictly regression.) Each neuron is connected to every neuron in adjacent layers but not to any others.

Following Lahav et al. 1996, each neuron $`j`$ in layer $`s`$ receives the $`N`$ outputs $`x_i^{(s1)}`$ from the previous layer $`s1`$ and gives a linear weighted sum over the outputs,

$$I_j^{(s)}=\underset{i=0}{\overset{N}{}}w_{ij}^{(s)}x_i^{(s1)}.$$
(9)
There is usually an additive constant, $`w_{0j}`$, where $`x_0=1`$, in this linear sum. This ‘bias’ allows the outputs to be shifted in analogy with a DC level.

The neuron then performs a non-linear operation (the transfer function) on the result to give its output $`x_j^{(s)}`$, typically a sigmoid or, as used here, the tanh function, which has an output range of $`1`$ to 1:

$$x_j^{(s)}=\frac{2}{1+\mathrm{exp}(2I_j^{(s)})}1.$$
(10)
The parameters are normalized to zero mean and unit variance. This is not strictly necessary as the net can in principle perform an arbitrary non-linear mapping, but it enables the weights to be initialized in the range $`1`$ to 1 and not be made unduly large or small relative to each other by the training. This is particularly helpful for larger networks.

The weights are prevented from growing too large by using weight decay, a regularisation method which adds a term $`d`$ to the cost function which penalizes large weights:

$$d=\mathrm{const}\times \frac{1}{2}\underset{j}{}w_j^2.$$
(11)
Regularisation is also helped by the normalization.

The weights are adjusted by the training algorithm. In galaxy classification this has typically been the well-known backpropagation algorithm (Werbos 1974; Parker 1985; Rumelhart, Hinton & Williams 1986) or the quasi-Newton algorithm (e.g. Bishop, 1995). The Matlab software allows the specification of which one to use from a number of choices including these. Here another algorithm popular in neural net research is used: the Levenberg-Marquardt method (Levenberg 1944; Marquardt 1963, also detailed in Bishop 1995). This has the advantage that it is very quick to converge to a minimum of the cost function, and it is able to cope with steep gradients in the parameter-cost function space by approximating gradient descent, and with shallow gradients by approximating Newton’s method. It is thought to be the fastest algorithm for networks of up to a few hundred weights and its implementation in Matlab further improves its performance.

Following the neural network toolbox documentation, the algorithm works by using the fact that when the cost function has the form of a sum of squares the computationally expensive Hessian matrix H can be approximated as:

$$\text{H}=\text{J}^T\text{J},$$
(12)
and the gradient is:

$$𝒈=\text{J}^T𝒆,$$
(13)
where J is the (much easier to compute) Jacobian containing the first derivatives of the network errors with respect to the weights and biases and $`𝒆`$ is a vector containing the network errors, where the network error is the network type minus the target type.

The algorithm then performs the update:

$$𝒘_{k+1}=𝒘_k[\text{J}^T\text{J}+\mu \text{I}]^1\text{J}^T𝒆,$$
(14)
where I is the identity matrix and $`\mu `$ is the ‘momentum’. A large $`\mu `$ approximates gradient descent and $`\mu =0`$ is Newton’s method. $`\mu `$ is given a large initial value so that gradient descent enables the area of the minimum to be found quickly. It is then decreased after each step where the cost function reduces, thus moving towards Newton’s method which is faster and more accurate near the minimum.

Matlab allows a number of adjustable parameters for the training. The default values were used. The parameters include:

epochs: maximum number of training iterations (100)

min\_grad: minimum gradient of the cost function ($`1.00\times 10^{10}`$)

mu: initial value of $`\mu `$ ($`1.00\times 10^3`$)

mu\_dec: amount to multiply $`\mu `$ by when the cost function is reduced by a step (0.1)

mu\_inc: similarly for when the cost function increases (10)

mu\_max: maximum $`\mu `$ value ($`1.00\times 10^{10}`$)

The criteria used for stopping training were epochs, min\_grad, and mu\_max, whichever was reached first. An explanation for mu\_max being used is that, whilst appearing indicative of a diverging solution, it is in fact showing that the algorithm is unable to make a further step to reduce the cost function. The algorithm only steps if a resulting reduction is found, so it tries progressively larger steps to search for this, until mu\_max is reached. One could also use as a stopping criterion a validation sample, in which the training is stopped if the cost function when the network at that stage of its training is run begins to increase. However, with Levenberg-Marquardt the minimum may be reached in very few iterations (e.g. ten or less), and with the large training samples used here the validation sample gives virtually the same value of the cost function as the training sample. There is little danger of overfitting because of the size of the training sample and the intrinsic spread in the galaxy properties. An exception may be a large network with the eyeball training sample (see §4.1).

In general the space of parameters and cost function may have arbitrarily many local minima. It is thus necessary to start with several random initializations of the weights (or ‘runs’) to avoid a poor local minimum giving spurious results. The results can then be viewed with the poorest networks down-weighted or ignored, or by using the median type. Here the median type is used because although very few runs will be significantly poorer than average, the ones that are may be by enough such that the mean is a worse measure than the median. The median type quoted in this paper is always taken from ten runs. The typical scatter between runs is found to be significantly less than the mean RMS spread of the network types about the targets.

The trained network is then applied to the test sample, and it is for this sample that the tabulated results are recorded. The training and test samples must be independent but the training sample must be representative of the test sample. Here the galaxies are given in a random order, the first half was used for training, and the second half for testing. For the eClass and redshift one eighth of the DR1 galaxies were used for training and 10,000 for testing. The samples had their outliers removed, and those with eClass and redshift targets were evened, using the methods described in §2, resulting in training and test samples of approximately 10,000 galaxies (8,501 and 9,801 for eClass; 10,132 and 9,801 for redshift). These samples are easily large enough to train and test the networks without using undue amounts of memory. The resulting eyeball samples of 674 (training) and 683 (testing) were not evened as this would make the samples too small using the method here. For eClass and redshift the network was simulated on the rest of the DR1 sample with outliers removed (79,769 galaxies for both targets).

Further details on neural nets can be found in Bishop (1995) and in the context of galaxy classification in Lahav et al. (1996).

## 4 Results

The networks were iterated over many parameter sets, architectures and random initializations of weights. The results are shown for the parameter sets for the network architectures 1 (single neuron) and 8:1 (8 neurons in a hidden layer) in Table 2. Some of the best sets (highest correlation between network output and target type/lowest root mean square difference between network output and target type; the one almost always corresponds with the other) were run on more architectures. These are shown in Tables 3 and 4. The architectures shown give reasonable execution times, since the Levenberg-Marquardt algorithm has memory requirements which scale as $`N^2`$ where $`N`$ is the number of weights. The largest number of weights used is in the hundreds.

### 4.1 Effect of Network Architecture

Tables 3 and 4 show that a network with a single hidden layer with a few neurons is adequate for the task of predicting these galaxy parameters using Sloan data. Thus many network runs could be used to get a good distribution of the assigned type for any particular galaxy. Beyond about ten hidden neurons there is little improvement and in fact the standard deviation of assigned types to individual galaxies from the multiple initializations, usually much less than the RMS between actual and target types, starts to increase. A network, e.g. hidden units of 8:1, is clearly better than a linear mapping, represented by a single neuron, and although in some cases the improvement in correlation/rms is not large the plot of network type versus target type (as in Fig.s 13) is a much smoother function of target type. The networks are almost certainly limited in their performance by intrinsic scatter in the training sample. This can be seen if the network is tested on the sample it has just been trained on – its performance is very similar. This also confirms the earlier statement that overfitting is unlikely with the 8:1 nets and sizes of training samples used. The increased spread in assigned types with larger networks may be indicative of overfitting, particularly with the eyeball type as the number of weights becomes comparable to the number of training examples. The results presented in the figures use the 8:1 architecture, for which this is not a problem.

### 4.2 Effect of Parameter Set

In general it seems that certain parameters are good for predicting the targets, but that if all the parameters are added in, the correlation improves over the subsets. The correlation is not improved by duplicating the few best parameters, so it would appear that genuine information is present in the less good parameters and it is adding these and not just increasing the size of the network which helps. We therefore use all parameters in generating the Figures. The model magnitudes used are those from the SDSS DR1 which have been found to be offset by up to 0.2 mag, but this does not matter here, since the training and test samples are affected in the same way. As expected, including the magnitudes as well as the colours adds little to the correlation as no significant new information is added.

### 4.3 Results for the Different Target Types

#### 4.3.1 Eyeball Morphological Type

Previous studies (Naim et al. 1995; Lahav et al. 1995) have shown that neural networks are able to reproduce human-assigned morphological classifications with the same degree of accuracy as another human expert, about 1.8 types in the $`5`$ to 11 T type range. Here the types are assigned in bins of 0.5 in the range 0 to 6. Fig. 1 shows the median network type versus target type for ten runs. The network gives correlations up to 0.93 with an RMS of 0.55, about 9 per cent of the range, or the same as the width of the bins for the types. Smoothing the training sample over the bins by adding random noise of half the bin width was also tried but this did not improve the correlation, as the bins are quite small relative to the range in targets.

#### 4.3.2 eClass Spectral Type

The ANNs are able to predict the eClass spectral type when trained on galaxies with spectra in the SDSS with a correlation of up to 0.95 and RMS of 0.06 (4 per cent) for the test sample in the range $`1`$ to 0.5. The results for the simulation on the rest of DR1 are shown in Fig. 2. The shape is not perfect – a plot of net type - target type versus target type is not precisely symmetrical about zero, but the sigmoid shape seen when the training sample is not evenly sampled (§2) is not as pronounced. The sigmoid shape has been seen previously, e.g. Naim et al. (1995), where the network ‘avoided the ends of the scale’.

#### 4.3.3 Redshift

The generality of the method means that any parameter can be trained on and predicted, hence a photometric redshift can be obtained (Fig. 3). The correlation and RMS are up to 0.93 and down to 0.02. The RMS is comparable to other photometric redshifts in the literature found using neural networks, e.g. Tagliaferri et al. (2002) and Firth, Lahav & Somerville (2003), and to those derived from SDSS data (Csabai et al., 2003).

## 5 Discussion

The main result is that the networks can predict morphological classifications, spectral types and redshifts of galaxies using just photometric parameters. This paper uses moderately sophisticated neural net techniques on a data set of unprecedented size and quality. There are many further techniques which could be used, and possibilities to try out. In particular there are many sophisticated ANN techniques which have been little used in astronomy but which may now be justified by the size of the datasets available. However, with the current data it is unlikely that they would make large improvements as the results are almost certainly limited by the intrinsic spread in the training samples, and one can never improve upon the training sample.

Possibilities include, at a basic level, varying the galaxy and neural net parameters used here, for example the number of random initializations or various Matlab parameters. More sophisticated neural net techniques include a better even sampling of the training sample so that targets are evenly spread over the range – here over-populated bins are simply cut down to size but under-populated bins are not altered. This may be especially useful for star formation rate and can be done using K-means clustering or a self-organising map (e.g. Tagliaferri et al., 2002). Improved regularisation, e.g. hierarchical Bayesian learning as opposed to weight decay, could be implemented. Multiple outputs for the network could be used to perform classification as opposed to regression. This is complementary work rather than an improvement and could be used for any of the types, in particular the eyeball types, or it could be used to e.g. assign probabilities to photometric redshift bins, as each output can give the a posteriori probability that the type is that output given its input parameters. This would also show objects for which the photometric redshift is less certain, as there may be no one bin with a high probability, or it may be split with peaks occurring in two separated bins. A different learning algorithm, e.g. quasi-Newton or conjugate gradient, would be needed for a classifier as Levenberg-Marquardt requires one output. Other learning algorithms could also be used with a validation sample. The disadvantage of a classifier is that the number of bins for the output is fixed for the network used. With the regression used here one can bin the assigned types afterwards if desired. Various methods exist for using committees of networks, apart from that here of using multiple random starts and using the median type assigned. Examples include constructive learning, bootstrap training samples, forward selection, backward elimination, cross validation and waterfall. There are other methods for global optimization apart from multiple random starts, e.g. simulated annealing and genetic algorithms. Many of these possibilities (regularisation, committees, etc.) are discussed at the comp.ai.neural-nets newsgroup FAQ at ftp://ftp.sas.com/pub/neural/index.html, whilst a waterfall of networks was used in galaxy classification by Adams & Woolley (1994). A further possibility is the use of unsupervised networks, i.e. those in which a predefined similarity criterion is used and the data is left to organize itself, with no training sample required. Unsupervised networks objectively find clusters of similar points in a dataset and can be used as a basis for classification. A well-known unsupervised network which has been used in classifying galaxies is the Kohonen self-organizing map (Kohonen, 2001), used by Naim et al. (1997b) for galaxy morphology. Naim et al. (1995) used principal component analysis to reduce a set of 24 galaxy parameters to 13 and found that the latter was as good for predicting types. This was tried, and found to be quicker for networks with more than 100 weights over ten runs but it was found that the correlations were generally slightly worse, as some information was lost (and it cannot be gained by PCA). The time taken was mainly that for the PCA, then the $`N^2`$ scaling for Levenberg-Marquardt with $`N`$ weights.

The usefulness of the methods here is that they are able to predict either spectral parameters using just photometry or assign morphological types at a vastly greater rate than humans but to the same accuracy. Much can be done with the types once they have been assigned, and this will form the basis of future work, in the distributions of these types and in their use to augment large scale structure studies using SDSS data with other physical measures such as colours.

The SDSS Southern Survey (York et al., 2000) is repeatedly imaging a smaller area of the southern galactic cap to go fainter in imaging and spectroscopy than the northern survey. Spectra from the Southern Survey could be used as training samples for galaxies at higher redshifts and below the northern spectroscopic flux limit.

One could thus look at galaxy evolution according to any assigned parameter. One could also, for example, project galaxies of unknown redshift about ones with known redshift, or push fainter down the luminosity function if assumptions are made about clustering. A particular statistic of interest is that of marked point processes (e.g. Beisbart, Kerscher & Mecke, 2002), in which the effects of intrinsic variation and those of environment can be separated.

Also, if one could predict physical parameters directly this would be extremely useful. One example is the star formation rate. The sample of 8,683 galaxies detailed in Gomez et al. (2003) was investigated. This is a volume limited sample from $`0.05z0.095`$ with well measured redshifts and H$`\alpha `$ star formation rate. At present the star formation rate is poorly predicted by the ANN, being best at zero but widely spread above this. Improved results may be obtained for networks trained on just those galaxies which are star forming.

Further targets which could be predicted include the bulge to disc ratio or the 2dF $`\eta `$ spectral type, and there are further galaxy parameters which could be used such as Sérsic indices, or spectral parameters for predicting morphological types.

It is not immediately obvious whether the resulting distributions say more about the galaxies or the assigned types, but with the numbers of galaxies available biases in the assigned types from the network could be studied in detail. Any biases of this sort are already less than the intrinsic spread in assigned type and one could compare results using a sample where the target types are available to see if different results are obtained. If not, then as long as the sample used has photometry of which the training sample was representative, the network types can be used with confidence.

## 6 Conclusions

The neural nets are able to predict the eyeball morphological type, the spectral type eClass, and the redshift using parameters available for all galaxy images in the Sloan Digital Sky Survey Data Release One. The correlations are 0.93, 0.95, and 0.93 respectively. The mean RMS errors between the network output and the known type for a set of unseen galaxies of which the training set formed a representative part are 0.55, 0.06 and 0.02 (approximately 9, 4, and 5 per cent of the ranges of the targets).

## Acknowledgments

Nick Ball thanks Andy Connolly, Andrew Hopkins, Ofer Lahav, Osamu Nakamura, Bob Nichol, Stephen Odewahn, Andrew Ptak, Kazuhiro Shimasaku, Alex Szalay, Chisato Yamauchi and the comp.ai.neural-nets and comp.soft-sys.matlab newsgroups for useful discussions and/or information.

Nick Ball is funded by a PPARC studentship.

Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/.

The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, University of Pittsburgh, Princeton University, the United States Naval Observatory, and the University of Washington.


text:
# Study of Equilibrium Using Collision Dynamics

## 1 Introduction

The intermediate energy heavy-ion collisions are very useful to study the non - equilibrium dynamics of finite size systems. In addition, one also has an opportunity to understand the properties of hot and dense nuclear matter that exists for a short span of time during the reaction. It is worth mentioning that no direct extraction can be made about these phenomena and therefore, one has to rely either on an indirect method of extraction or on a dynamical theoretical model that is capable of simulating the reaction from the start till the end where measurements are made. The question of equilibrium (and thermalisation) can be dealt within those models which are not based on the assumption of any equilibrium. The transport models such as Boltzmann-Uehling-Uhlenbeck \[BUU\]  or Quantum Molecular Dynamics \[QMD\]  are very helpful as these models can also handle the non-equilibrated nuclear matter formed during the early phase of the reaction. Naturally, this anisotropy in the momentum space (that reflects the non-equilibrium situation) should be taken into account while studying a heavy-ion collision.

We here plan to discuss the degree of equilibrium that can be reached in an intermediate energy heavy ion collision. This study is made within the framework of QMD model . We shall show that the degree of equilibrium (studied via rapidity distribution) depends on the reaction geometry as well as on the number of collisions any individual nucleon suffers. The preservation of the initial memory of nucleons is directly related with the number of collisions one suffers. Therefore, the momentum space of those nucleons who suffer large number of collisions should be better thermalized.

The section 2 deals with details of the model. Section 3 depicts the results and discussion. Finally, we summarise the results in section 4.

## 2 The Model

The nucleons in a molecular dynamics picture interact via two- and three-body forces. The explicit two- and three- body interactions lead to the preservation of fluctuations and correlations that are important for N-body phenomena like multifragmentation. This is in contrast to the one- body dynamical models which are suitable for one-body observable only.

In QMD model , each nucleon is represented by a Gaussian distribution whose centriod propagates with classical equations of motion:

$$\frac{d𝐫_i}{dt}=\frac{dH}{d𝐩_i};$$
(1)
$$\frac{d𝐩_i}{dt}=\frac{dH}{d𝐫_i},$$
(2)
where the Hamiltonian is given by :

$$H=\underset{i}{}\frac{𝐩_i^2}{2m_i}+V^{tot},$$
(3)
with

$$V^{tot}=V^{loc}+V^{Yuk}+V^{Coul}+V^{MDI}$$
(4)
Here $`V^{loc},V^{Yuk},V^{Coul}`$ and $`V^{MDI}`$, represent, respectively, the Skyrme, Yukawa, Coulomb and momentum dependent (MDI) parts of the interaction. The interaction without MDI part is called static interaction. The different values of the compressibility in the Skyrme force give possibility to look for the role of different equations of state termed as soft and hard equations of state. The inclusion of momentum dependent interactions are labelled as soft momentum dependent (SMD) and hard momentum dependent (HMD), respectively.

The G- matrix at higher excitation energies becomes complex in nature and its imaginary part acts like the collision term. We shall use here an energy dependent nucleon-nucleon cross-section. It is, however, worth mentioning that the reaction dynamics depends on the form and magnitude of the nucleon-nucleon cross-section .

## 3 Results and Discussion

The present study is made by simulating the reactions of Ca-Ca, Xe-Sn and Au-Au at different incident energies as well at different impact parameters. We here use a soft equation of state along with energy dependent nucleon-nucleon cross-section  through out the discussion.

There are several different ways to define the degree of equilibrium. The first quantity is the anisotropy ratio $`R_a`$ which is defined as

$$R_a=\frac{\sqrt{p_x^2}+\sqrt{p_y^2}}{2\sqrt{p_z^2}}.$$
(5)
This anisotropy ratio is an indicator of the global equilibrium of the system. The word global is due to the fact that this quantity does not depend on the local positions of nucleons and therefore, represents the equilibrium of the whole system. A full global equilibrium demands the anisotropy ratio to be close to unity. Another way to study the local equilibrium is to look for the relative momentum of two colliding Fermi spheres which indicates the deviation from a Fermi sphere. Note that the concept of local equilibrium is used by the hydrodynamical models to simulate the heavy-ion reactions. We shall, however, address the question of equilibrium with the help of rapidity distribution that also shows the stopping of nuclear matter in heavy-ion collisions.

The rapidity distribution can be defined as :

$$Y(i)=\frac{1}{2}\mathrm{ln}\frac{𝐄(i)𝐩_z(i)}{𝐄(i)𝐩_z(i)},$$
(6)
Where $`𝐄(i)`$ and $`𝐩_𝐳(i)`$ are, respectively, the energy and the longitudinal momentum of the i-th particle. For a full equilibrium, one should get a Gaussian shape distribution peaked at the mid-rapidity region. In other words, both the anisotropy ratio and rapidity distribution are related and can give insight into the equilibrium process of a reaction.

Using the above description, we plan to relate the frequency of nucleon-nucleon collisions with the degree of equilibrium. This is achieved by dividing the rapidity distribution of the whole system into different collision zones. Among such zones, the spectator matter (SM) consists of all those nucleons who do not suffered any collision. We divide the participant region into three zones : (i) the low collision matter (LCM) which contains all nucleons who suffer less than four collisions. (ii) the moderate collision matter (MCM) that has those nucleons with 5-9 collisions and (iii) high collision matter (HCM) which takes care of all those nucleons with more than ten collisions .

First of all, we study the evolution of rapidity distribution. In figure 1, we display the evolution of rapidity distribution dN/dY for the reaction of Xe-Sn at 400 MeV/nucleon for different geometry. Here we display the evolution for central (b=0 fm), semi-central (b=4 fm) and peripheral (b=8 fm) impact parameters.

The central reaction leads to a very high nucleon-nucleon collision rate and density whereas the peripheral reaction has a very small overlap. At the beginning, whole matter is grouped into either projectile matter or into target matter. No nucleon-nucleon collisions occur till 9 fm/c , therefore, no change in the rapidity structure appears. Due to some nucleon-nucleon collision between 9 and 18 fm/c, few nucleons shifts to the mid-rapidity zone. On the other hand, drastic changes can be seen (for the central and semi-central reactions) between 18 and 28 fm/c. It is that time when the density and collision rate (not shown here) is maximal. The low frequency of the nucleon-nucleon collisions in peripheral geometry does not allow substantial changes in the shape of the rapidity distribution even at the end of the reaction. This also points towards the lack of degree of equilibrium in peripheral collisions. On the other hand, a nearly equilibrium can be seen in central collisions suggesting better equilibrated nuclear matter in these reactions. The study of anisotropy ratio (not shown here) also depicts similar picture. Its final values are close to 0.87, 0.67 and 0.37, respectively, for central, semicentral and peripheral reactions. As mentioned in the beginning, for a complete equilibrium, this ratio should be equal to one. Both the rapidity distribution and the anisotropy ratio suggest that central reactions should be better thermalised. It is worth mentioning that both these quantities depends on the mass of the system .

It will be of further interest to analyse the momentum space of nucleons of equilibrated matter. This is achieved by dividing the final state rapidity distribution into the above mentioned collision zones (see figure 2).

Naturally, the frequency of the nucleon- nucleon collisions increases with the incident energy, therefore, HCM percentage is much more at 1 GeV/nucleon than at 400 MeV/nucleon. The maximal collision rate for Xe-Sn reaction at 400 MeV/nucleon is 34 whereas it is 56 at 1 GeV/nucleon. We find that the gaussian shape is more prominent for those nucleons that suffer at least ten collisions. This is true at both incident energies. The nucleons with less than 10 collisions observe a partial equilibrium.

In other words, a complete equilibrium is possible only for those nucleons that suffer at least 10 collisions in a reaction. The nucleon- nucleon collisions destroy the initial correlations, therefore, are dominant mode of achieving the thermalization in a reaction. From figure 2, it is also evident that the momentum distribution of nucleons with less than 4 collisions is still like that of projectile and target.

The mass dependence in heavy ion reactions has been found to affect the dynamics ranging from the fusion (at low incident energy) to the multifragmentation as well as collective flow and particle production (at intermediate and relativistsic energy) . To see the role of mass dependence in thermalization, we display in figure 3, the rapidity distribution for different collision zones using the semi-central reactions of Ca-Ca, Xe-Sn and Au-Au at 400 MeV/nucleon. The conclusions drawn from figure 2 are also valid here. We find that independent of the mass of the system, at least ten collisions are needed to form nearly equilibrated nuclear matter. It is also evident that the heavier matter are better thermalized compared to lighter systems. Similar trends can also be seen from the evolution of the anisotropy ratio. In other words, the degree of equilibrium depends on the size of the interacting system.

In figure 4, we display the transverse momentum $`(Px/A)`$ as a function of rapidity distribution for the semi-central reactions of Ca-Ca, Xe-Sn and Au-Au at 400 MeV/nucleon. We see that the spectator matter does not have any transverse momentum whereas the intensity of transverse flow increases with collision rate. The HCM has maximum transverse flow. Further the transverse flow depends strongly on the mass of the system which is in agreement with all previous observations and calculations . If one extrapolates the above results to very low incident energy, one will find that the balance energy ( i.e. the energy at which flow disappears) will be lower for heavier systems compared to lighter systems.

## 4 Conclusions

Summarizing, within a dynamical quantum molecular dynamics model, we present the study of equilibrium (i.e. thermalization) in heavy ion collisions at intermediate energies. This was achieved by simulating the reactions of Ca-Ca , Xe-Sn and Au-Au at different geometry as well as incident energies. We find that the momentum space of nucleons that suffer at least ten collisions is better thermalization whereas nucleons with less number of collisions exhibit a partial equilibrium. The nucleonic flow also shows strong collision dependence.

## Acknowledgement(s)

This work is supported by the Department of Science and Tehnology, Government of India.