text:
# Securing a HENP Computing Facility

## I Background

The RHIC Computing Facility (RCF) at the Brookhaven National Laboratory (BNL) is the primary data archiving and processing facility for the Relativistic Heavy Ion Collider (RHIC) at BNL. The US ATLAS Computing Facility at BNL is the US Tier 1 computing facility for the ATLAS experiment at the Large Hadron Collider (LHC) at CERN. As major computing centers, these facilities provide a myriad of services to users that are on-site at BNL and off-site at collaborating institutions. These services include the archive storage of raw and analyzed data from the experiments, interactive and batch data processing services and general interactive login facilities. Traditionally, such facilities have been completely open to the Internet and have provided services, e.g. FTP (File Transfer Protocol), telnet, IMAP (Internet Message Access Protocol), that were not designed with cybersecurity in mind.

An intrusion into either facility by a malicious intruder can result in damages from which recovery can be extremely costly in time and money. Some of the difficulty is associated with the enormous facility “state” represented by the 100’s of terabytes of on line disk storage and the 1000’s of compute nodes at the facility. However, this is not the limits of the damage. Since the facility is part of the data collection process, disruption of the facility can also result in the loss of data taking ability at each experiment. This loss is extremely costly since operating an accelerator is an expensive undertaking. Finally, because the facility is the sole archive of a substantial fraction of raw and analyzed data created by the experiments, loss of these archives from actions by malicious intruders can result in irreparable damage to on going scientific research.

Although many other computing facilities are at risk from cybersecurity incidents, two properties make the RCF and US ATLAS computing facilities high profile targets. First, the vast amounts of computing resources make the facilities appealing targets to intruders. Second, the location of the facility within a US government, Department of Energy facility, with its associated .gov top level domain makes these centers choice targets for politically motivate malicious intruders.

With the proliferation of Internet connectivity outside of the academic and research communities, the open data center architecture needs to be modified to mitigate the increased risk of cybersecurity threats.

## II Securing the Facility

The process of improving the security of the facilities was a multi-stage process. The first step involved the identification of assets and services at the facilities. In the second step, priorities were assigned to the identified asset or service with respect to importance. The third step was identifying the potential threats to the assets or services. The fourth step was modifying the network topology of the facility to mitigate the potential threats. The fifth step was replacing insecure services with more secure services. The sixth step was the modification of user and administrator behavior to increase the security of the facility. The seventh step was an assessment of the effectiveness of the instituted changes. The final step was the iteration of the proceeding seven steps to continue to improve the security of the facility.

### II.1 Assets and Services

There are a number of assets and services at the RCF and US ATLAS computing facilities with differing security needs. The main data store, the High Performance Storage System (HPSS) hpss-ref , contains the raw and analyzed data that is generated from the experiments. The protection of the stored data and the operation of the service are paramount to the facility. Without it, the other components of the facility are inconsequential. The service provide by HPSS is data storage and retrieval via an FTP-like mechanism. The operation of the system is effectively independent of the other assets and services in the facility.

After HPSS, the next service of importance at the facility is the authentication and authorization infrastructure, consisting of a cluster of Network Information System (NIS) servers. Without account information, the compute farm and the NFS (Network File System) storage are unusable. As a result, the protection of this information is crucial.

At approximately equivalent levels of importance are the NFS servers and the Linux based computation farm. The NFS storage consists of approximately 150TB of disk space and 20 NFS servers. Although, extremely disruptive, loss of NFS data or service is not catastrophic. Lost unprocessed data can be restored from HPSS and lost processed data can be recovered by reprocessing the raw data. However, depending on the details of the data loss, the recovery process may be short and simple or long and complex. The Linux based computation farm consists of approximately 1100 dual CPU systems and 130TB of local storage. Loss of these systems would halt data processing and loss of data stored on the local disks would require reloading of data from the data source or the recreation of data by reprocessing. As with the NFS servers, the type of loss would dictate the cost of recovery.

The wide area accessible Andrew File System (AFS) openafs-ref  storage is at the next level of importance at the facility. Loss of this service would be inconvenient for data processing on site. Recovery of data corrupted AFS storage would entail restoration from backups, resulting in a maximum loss of about 24 hours of work. Loss of AFS service would definitely cause major disruption to work at off site locations. For on site users, loss of AFS service can, if necessary, be supplanted with NFS service.

The remaining services provided by the RCF and US ATLAS facilities include web service, Samba service, and email service. Although important for the users, these services are not coupled to the higher priority mission of the facilities, which is to provide storage and computational resources to the individual experiments.

With the assets and services at the facilities identified, an assessment of the security threats is necessary to determine the changes that need to be made to increase the security of the facilities.

### II.2 Assessing the Threats

The threats to computing resources at the RCF and US ATLAS computing facilities fall into two classes, direct network assaults and compromised accounts. Examples of direct network assaults include direct attacks on network service daemons (e.g., buffer overflow exploits, protocol weaknesses), hijacking of network connections (e.g., replay attacks, forged packets), denial of service (DoS), distributed DoS, and web based attacks (both client and server). This class of attacks do not require login access on a facility system. Examples of compromised accounts include stolen or cracked passwords, social engineered access, dead accounts, and malicious insiders. This class of attacks involve gaining access to a login account at the facility from which exploits are launched. An interesting distinction between these two classes of assaults is the mitigation of network attacks can be more easily accomplished with technology and architecture changes. On the other hand, behavioral changes on the part of the system administrators or users are usually required to fix attacks in the compromised accounts class.

### II.3 Facility Architecture

A modification of the facility architecture is the simpler of the two classes of changes that were made at the RCF and US ATLAS facilities. (The harder change being the modification of user and administrator behavior.) The two main goals of the architectural modifications were compartmentalization and isolation. Where possible, unrelated assets and services were separated so that security breaches in one location did result in immediate security breaches in another location. Also, where possible, network accessability to assets and services was minimized.

The core facility, consisting of the HPSS servers, NFS servers, AFS servers, NIS servers and Linux computation farm, were placed behind a firewall with a “default deny” firewall rule. With the exception of the AFS servers and AFS service, there is virtually no direct network access from the outside to the core facility. In this configuration, the experiment counting houses, where raw experiment data is generated, are placed within the firewall.

For scalability and availability, the facility firewall was a multi-component firewall, consisting of a Cisco PIX firewall, augmented by ssh openssh-ref , secure file transfer and Samba samba-ref  (SMB) gateways. The gateways were implemented using general purpose computers running host based firewalls (ip-filter ip-filter-ref  or iptables iptables-ref ) on stripped down and hardened operating system installations. Each gateway system is equipped with two network interfaces (NIC), allowing them to bypass the PIX firewall, thus providing some scalability and availability. This dual NIC configuration is what dictated the need for the host based firewalls. The separation of the ssh and data transfer gateways into separate systems simplifies the hardening and maintenance of each system. It also isolates interactive login traffic with its low latency requirements from data transfers which are relatively insensitive to latency and typically result in high system load.

In conjunction with the introduction of the gateways was the elimination of all protocols with clear text passwords. FTP data transfers were replaced with scp openssh-ref , bbftp bbftp-ref , and ssh protected ftp. (The latter being normal FTP with the FTP control channel being sent through an encrypted ssh tunnel.) Interactive logins via telnet and rsh were replaced with encrypted logins via ssh. This latter change providing the additional benefit of protecting X x-ref  window traffic.

With the new core facility configuration, triage scenarios are possible with respect to intrusions. The first level would involve disconnecting the facility from the Internet, thus allowing on site user to work and on going operations to continue. The second level involves disconnection from the on site network. This would allow on going operations to continue, including data taking, but would cut off on site and off site users. The third level would involve the shutdown of all servers and services except for HPSS, thus halting all operation except data taking.

Located both outside and inside the facility firewall are multiple web servers. The services provided and resources used by the web servers are carefully distributed to the appropriate web server. Located outside the firewall is a non authenticating web server, serving both static and dynamic (aka CGI scripts) pages. This web server is completely stand alone, requiring no services or resources from systems inside the firewall. Two web servers are located inside the firewall, an authenticating web server and a non authenticating web server. The authenticating web server serves both static and dynamic web pages, but requires users to authenticate to the server before access is granted to the pages. (The authentication occurs over an SSL protected connection openssl-ref .) The assumption is that authenticated users, i.e. legitimate facility users, are not malicious and will not attempt to exploit weakness that may be present in the dynamic web pages. This web server is located behind the facility firewall because of its dependence on a core facility service. The second web server located inside the firewall is a non authenticating web server that serves only static user pages. Since it depends on access to user home directories, it resides behind the firewall. To eliminate the problems associated with CGI scripts and other dynamic content, only static pages are served.

The final component of the new facility architecture was the separation of email services from the core facility. A standalone SMTP and secure IMAP server imap-ref  was installed outside of the facility firewall. The email server utilizes a separate password database that is independent of the password database used for interactive access to the core facility. This separate database can potentially isolate the core facility from email account compromises. In addition, the use of SSL protected IMAP (replacing standard IMAP) protects email passwords while they traverse the network.

### II.4 Human Factors

Architectural changes to the RCF and US ATLAS facilities combined with technology choices significantly enhances the security of the facilities. However, a substantial amount of the security is dependent on the actions and behavior of users and administrators. As an example, benefits of encrypted ssh connections to the facility are negated if users explicitly set the X window DISPLAY variable to bypass the encrypted tunnel, or use telnet to log into the system from which he/she runs the ssh client. Similarly, system administrators can destroy site security if they are “social engineered” into providing account access to unauthenticated or unauthorized individuals.

Change user and administrator behavior primarily consists of education to make people aware of the security issues and to provide information on secure ways of doing tasks. However, education is not exceedingly effective. There is a distinct difference between awareness and understanding and theory and practice. Some behavior was enforced, for example disabling of telnet and proactive password checking with npasswd npasswd-ref . Unfortunately, most behavior cannot be enforced.

## III Operational Assessment

Operational experience that has been collected with the new facility architecture in place provides information about what is working and what is not. As was expected, migration to the new configuration was a long and drawn out process, much akin to changing the direction of an oil supertanker. New tools and configurations needed to be tested, documented and put into production. In addition, generous lead time was required so that users could prepare for the system changes. At this point, the facility has accomplished the first iteration of the facility hardening process outlined at the beginning of the previous section.

The critical components of the new facility are in the multicomponent firewall, since they represent the interface to the facility for the users. The most frequently used component in the firewall are the ssh gateways. Implemented using four 650MHz Intel Celeron PCs, they are each able to handle up to 80 interactive users running xterms. However, a single graphics intensive application can saturate the CPU, resulting in longer latencies for users. Additional problems have been encountered with scp data transfers on the ssh gateways. As with graphics intensive applications, a single scp session can saturate the CPU resulting in responsiveness issues for other users. The ssh gateways have subsequently been upgraded to 1.7GHz Pentium 4 systems, providing substantially improved performance.

The FTP gateway, currently a dual 800MHz Pentium III system with dual gigabit ethernet interfaces, has been a little problematic. Getting users to use the FTP gateway was the first problem that was encountered. Over time, this has become less of an issue. An on going problem is the use of sftp and scp for bulk data transfers. As on the ssh gateways, a single sftp/scp transfer can saturate a single CPU. For a handful of power users, bbftp has worked well for bulk data transfers since only the control channel is encrypted. However, its non FTP-like user interface and its relative obscurity compared to scp and FTP are a barrier to wide spread adoption by users. Additional problems are associated with the use of NFS mounted disks as a source and destination for data files on the gateway. The performance of the NFS disk is a bottleneck for high speed transfers of individual files. The installation of local disk is an option but creates problems with disk space management. To alleviate some of the performance issues, the FTP gateway will be replaced with two dual 2.4GHz Pentium 4 systems.

For all components of the firewall, configuration and maintenance is problematic. The maintenance of the rules in the PIX firewall to handle the addition, removal, and movement of services is a significant undertaking. Judging the security of new service is particularly problematic. In addition, some services are difficult to protect with a firewall. With the host-based firewalls, the configuration of the rules can be tricky and can result in unexpected problems.

Finally, significant work is still necessary on audit trails and intrusion detection on the various firewall components.

## IV Future Work

Operational experience has revealed areas where modifications can be made to improve site security. In addition, the Internet and cybersecurity threats are constantly changing, as a result, defenses need to be upgraded to handle them. Since most of the security features in the new facility architectures are designed to handle the threats circa the late 1990’s, work needs to be done to deal with new types of threats.

One area of work is the additional hardening of the “edge” systems, i.e., those systems that are directly accessible from the Internet. Better understanding of the issues and technologies makes it possible to better protect the systems and the facility. Examples of these changes include better firewall rules, updated software, and additional re-architecting of services to improve security and manageability.

With the proliferation of sophisticated web services, the protection of web clients from malicious web servers is also an idea to be considered. On a similar vein, proactive institution of firewall rules to mitigate the risk of the facility being a launch site of a network attack is also worth considering.

### IV.1 Grid Services

Looming on the horizon are issues presented by Grid globus-ref  technology. With the Grid comes new types of security concerns. One promise of Grid technologies is the availability of global computation resources from virtually anywhere in the world. However, this also adds a new scale to the ramification of security breaches, e.g., access to global computational resources to launch distributed DoS attacks and the destruction or disabling of resources on a global scale. With the Grid, responsibilities, facilities and authorities become geographically and organizationally distributed. Command, control and communication become much more difficult.

Operationally, the phasing in of incomplete or untested Grid services into a production facility has unknown risks. Additionally, the skill sets and experience need to manage these new services do not exist and will likely result in operational errors which may have dire security consequences.

## V Conclusions

With a modification of the architecture of the RCF and US Atlas Tier 1 computing facilities, the replacement of insecure services with more secure services, and the education of both facility users and administrators, substantial gains in facility security have been obtained. However, security is an on going process, not an endpoint. The introduction of new requirements and services are making security harder and security breaches potentially more disastrous. Maintaining the integrity of the current security profile while growing the facility and accommodating new requirements and services will make security an increasingly difficult task.

###### Acknowledgements.

The authors wish to thank the five RHIC experiments, the US ATLAS group, the BNL Physics department, and the Network group in the Information Technology Division at BNL.


format: text

text:
# Deuterium fractionation on interstellar grains studied with the direct master equation approach

## 1 Introduction

Rate equations have long been used in modelling the surface chemistry that occurs on cold dust particles in the interstellar medium (Pickles & Williams, 1977; Hasegawa, Herbst & Leung, 1992). It has been pointed out (Charnley, Tielens & Rodgers, 1997) that, in the case when the average number of reactive species on a grain is sufficiently small, the rate equations may not be suitable to describe the system, since they do not take into consideration its discrete nature. One alternative is to use a Monte Carlo approach (Tielens & Hagen, 1982; Charnley et al., 1997; Caselli et al., 2002; Charnley, 2001).

Recently another approach to this problem (Green et al., 2001; Biham et al., 2001), has been proposed. Like the Monte Carlo method, this approach is based on a master equation. But, unlike the Monte Carlo technique, the master equation is directly converted into differential equations. In particular, the rate equations for the population of highly reactive species with small surface abundances are replaced by sets of differential equations for the probability that certain number of these atoms or molecules exist on the surface at the same time. One advantage of this new exact approach is its easy coupling with the differential rate equations used to model the gas-phase chemistry. A comparison of the two approaches is to be found in Stantcheva, Shematovich & Herbst (2002).

We have, so far, successfully applied the direct master equation approach to a system consisting of H and O atoms and CO molecules accreting onto a grain surface, and reacting to form a variety of reaction products (Stantcheva et al., 2002). In this paper, we report an expansion of the system to include D atoms accreting onto the surface and reacting with the other reactive surface species to produce deuterated isotopomers. The system is identical to that studied in Caselli et al. (2002) with modified rate equations and a Monte Carlo approach. A somewhat smaller network of reactions was utilised by Charnley et al. (1997). As opposed to Caselli et al. (2002), we concentrate on the doubly, triply and quadruply-deuterated species which we believe are likely candidates for observation in highly fractionated low-mass protostellar sources (Parise et al., 2002). The main purpose of this paper; however, is to show how the direct master equation method can be used efficiently.

The paper is organised as follows. First, in Section 2 a short description of the system under investigation is given. The theory, including the formulae used in the calculations, is given in Section 3. Section 4 presents the results, and is followed by a discussion in Section 5.

## 2 The chemical network and physical conditions

We have considered a system in which H, O, and D atoms and CO molecules accrete onto the grain surface and react to form the stable molecules H2, HD, D2, H2O, HDO, D2O, H2CO, HDCO, D2CO, CH3OH, CH3OD, CH2DOH, CH2DOD, CHD2OH, CHD2OD, CD3OH, CD3OD, O2, CO2, and the highly reactive radicals OH, OD, HCO, DCO, H3CO, H2DCO, HD2CO, and D3CO. The system pertains to old dense interstellar clouds, in which hydrogen is mostly molecular and carbon is mostly in the form of CO. The calculations were performed mainly at a temperature T=10 K, for a period of 104 years, and fixed gas-phase abundances of the accreting species, an approximation justified by the relatively short period of evolution of the surface chemistry compared with the time scale of the gas-phase chemical processes. At times even earlier than 104 yr, a steady-state condition is reached in which the surface populations of the stable species grow linearly with increasing time while the reactive species have fixed populations. The number of available sites per grain, which helps to define the diffusion rates, was taken to be $`10^6`$; this number refers to so-called classical grains with a size of 0.1 $`\mu `$m.

Three different sets of values for the gas-phase densities of H, O, CO, were used throughout the calculations. Given in Table 1, these sets of values are referred to as the low, intermediate, and high density cases because they were derived from gas-phase models at total densities of 103, 104, and 105 cm-3, respectively. In all cases, the concentration of H is near 1 cm-3. Unless stated to the contrary, all results are for a gas-phase abundance for D of 0.3 cm-3. This very high value relative to the atomic H abundance is presumably produced via fractionation in the cold gas, although current gas-phase models cannot easily reproduce such a value (Roberts, Herbst & Millar, 2002). Variations in the atomic deuterium abundance are also considered.

In calculating the accretion, desorption, and diffusion rates of the surface species, we followed the methods used in Hasegawa et al. (1992) – see also Caselli et al. (2002, 1998) – and used their values for the parameters necessary for the calculations. All particles were considered to diffuse over the surface solely via thermal hopping except for H and D atoms, for which quantum tunnelling was also considered. The rates used for diffusion are the so-called fast rates (Stantcheva et al., 2002) because these magnify the differences between exact and approximate methods for studying diffusive surface chemistry. Table 2 shows the accretion, evaporation, and diffusion rates over an entire grain for the accreting species in the model. Note that at the temperatures considered, the other heavy species in the model do not diffuse or evaporate at non-negligible rates.

The surface reactions used in our model as well as any non-zero activation energies $`E_\mathrm{a}`$ are given in Table 3. More details can be found in Caselli et al. (2002).

## 3 Master equation for the system

In the direct master equation method, the system corresponding to the grain surface is represented by a multitude of states that represent different possible populations of the surface species. In one state we might have 0 particles of the species A, 1 particle of the species B, 2 particles of the species C, …, whereas a different state might consist of 1 particle of A, 0 of B, 2 of C, …. With each state, we associate a probability for the system to be in this state, and we develop equations for the rate at which these probabilities change (the master equation). Basic to the approach is the fact that the populations of some of the surface species are correlated. As a consequence, the method, at least formally, requires the consideration of the evolution of the system as a whole, as opposed to the rate equation approach, in which the evolution of the average population of any given species is followed separately.

In the most general case, all the surface species are taken into consideration in each state of the system. Since the standard rate equations, however, give a sufficiently accurate description of the evolution of high-abundance species (Stantcheva et al., 2002), the inclusion of these species in the master equation will impose an unnecessarily heavy load on the computing resources given the need to integrate many coupled equations simultaneously. Thus, a more practical approach involving fewer coupled equations is to include only the highly reactive and low-abundance surface species in the master equation, and to use rate-like equations to solve for the surface populations of the rest of the species. As a guide to determine which species should be included in the master equation, one can use the results obtained via rate equations and, if available, Monte Carlo simulations. More generally, these species are reactive atoms and radicals, especially those that are precursors of major grain species. Here we include the 11 species H, O, D, OH, OD, HCO, DCO, CH3O, CH2DO, CHD2O, and CD3O. We enumerate these from 1 to 11 and refer to them as “probabilistic”. The remainder of the species are referred to as normal and are treated via rate-like equations.

Once we have determined the probabilistic species, we assign a probability $`P(i_1,\mathrm{},i_{11})`$ to any state {$`i_1,\mathrm{},i_{11}`$} that consists of $`i_1`$ H atoms, $`i_2`$ O atoms, $`i_3`$ D atoms, $`i_4`$ OH molecules, …, $`i_{11}`$ CD3O molecules, and solve the master equation, which consists of equations for the rate of change of the state probabilities $`P(i_1,\mathrm{},i_{11})`$. These equations are of the form

$$\begin{array}{c}\frac{dP}{dt}(i_1,\mathrm{},i_{11})=\hfill \\ \\ \underset{\{\mathrm{X}\}}{}k_{\mathrm{acc}}(\mathrm{X})n(\mathrm{X})\left[P(\mathrm{},i_j1,\mathrm{})P(\mathrm{},i_j,\mathrm{})\right]\hfill \\ \\ +\underset{\{\mathrm{X}\}}{}t_{\mathrm{evap}}^1(\mathrm{X})\left[(i_j+1)P(\mathrm{},i_j+1,\mathrm{})i_jP(\mathrm{},i_j,\mathrm{})\right]\hfill \\ \\ +\underset{\{\mathrm{X},\mathrm{Y}\}}{}k_{\mathrm{X},\mathrm{Y}}(i_j+1)(i_k+1)P(\mathrm{},i_j+1,\mathrm{},i_k+1,\mathrm{})\hfill \\ \\ \underset{\{\mathrm{X},\mathrm{Y}\}}{}k_{\mathrm{X},\mathrm{Y}}(i_j)(i_k)P(\mathrm{},i_j,\mathrm{},i_k,\mathrm{})\hfill \\ \\ +\underset{\{\mathrm{X}\}}{}k_{\mathrm{X},\mathrm{X}}\frac{(i_j+2)(i_j+1)}{2}P(\mathrm{},i_j+2,\mathrm{})\hfill \\ \\ \underset{\{\mathrm{X}\}}{}k_{\mathrm{X},\mathrm{X}}\frac{i_j(i_j1)}{2}P(\mathrm{},i_j,\mathrm{}).\hfill \end{array}$$
(1)
In eqs. (1), $`i_j`$ corresponds to the surface abundance of the species X, $`i_k`$ to that of the species Y, $`n`$(X) is the gas-phase abundance of X, and $`k_{\mathrm{acc}}`$(X), $`t_{\mathrm{evap}}^1`$(X), and $`k_{\mathrm{X},\mathrm{Y}}`$ are the accretion and evaporation rates of X, and the rate coefficient for reaction between X and Y, respectively (Hasegawa et al., 1992). In the units used here (Caselli et al., 1998), the rate coefficient is simply the sum of the the diffusion rates of X and Y over the entire grain. The first two terms on the right-hand side of eq. (1) represent the changes in the probability with time due to accretion and evaporation processes. The subsequent terms account for the changes due to reactions between two probabilistic species, both non-identical and identical. It is important to note that there should be two more terms on the right-hand side of eq. (1), which are present in our calculations but are not given here for considerations of simplicity. The first of these terms takes into account the change of the probability due to reactions between a probabilistic and a normal species, such as H+CO$``$HCO. The second accounts for (slow) reactions between two normal species, which have at least one probabilistic species as a product. In our current model there are no reactions of the latter type. In addition to these simplifications, we have also not indicated changes to more than two stochastic species in any term although there are reactions in which the populations of three such species can change (e.g. O + HCO).

The solution of eqs. (1) gives the probabilities for all the states, which, in turn, allows the calculation of the average abundance, or number of species per grain, of each probabilistic species, $`N_\mathrm{X}`$, as well as any necessary correlation terms, of the form $`N_\mathrm{X}N_\mathrm{Y}`$. These entities, then, are substituted into the equations used to solve for the abundances of the normal species. Although the latter equations are very similar to the widely used rate equations, they bear important differences. Equations (2) and (3) are the expressions for the average abundances of H2 and H2CO (both being normal species):

$`{\displaystyle \frac{dN_{\mathrm{H}_2}}{dt}}`$ $`=`$ $`t_{\mathrm{evap}}^1(\mathrm{H}_2)\times N_{\mathrm{H}_2}`$ (2)
$`+\mathrm{\hspace{0.33em}0.5}k_{\mathrm{H},\mathrm{H}}\times N_\mathrm{H}(N_\mathrm{H}1),`$
$`{\displaystyle \frac{dN_{\mathrm{H}_2\mathrm{CO}}}{dt}}`$ $`=`$ $`+k_{\mathrm{H},\mathrm{HCO}}\times N_\mathrm{H}N_{\mathrm{HCO}}`$ (3)
$`k_{\mathrm{H},\mathrm{H}_2\mathrm{CO}}\times N_\mathrm{H}\times N_{\mathrm{H}_2\mathrm{CO}}.`$
It can be seen that, unlike regular rate equations, both correlated abundances and terms in which 1 is subtracted from an abundance appear. The rate-like equations are obtained by summation over single and correlated probabilities for the relevant stochastic species involved in the formation and depletion of major species (Biham et al., 2001; Green et al., 2001; Stantcheva et al., 2002).

To propagate the system forward in time, both the master equation and the rate-like equations must be solved simultaneously. Since the probabilistic species have very low abundance, the probability for the system to be in states that contain high abundances of these species must be very small. It is reasonable, then, to neglect such states and limit our consideration to those that comprise only very few particles. For this purpose, we choose a set of parameters, {$`N_1`$,$`N_2`$,…, $`N_{11}`$} and $`N_{\mathrm{tot}}`$, such that we neglect every state, {$`i_1`$,$`i_2`$,…,$`i_{11}`$}, for which the following conditions hold true:

$`i_j`$ $`>`$ $`N_j,\mathrm{for}\mathrm{any}j,\mathrm{or}`$ (4)
$`{\displaystyle \underset{j}{}}i_j`$ $`>`$ $`N_{\mathrm{tot}}.`$ (5)
Because the first three particles in {$`i_1`$,$`i_2`$,…,$`i_{11}`$} correspond to H, O, and D, we require that $`N_1`$,$`N_2`$, and $`N_3`$ be equal to at least 2 so that molecular hydrogen, oxygen and deuterium can be produced. Thus, the set of minimum values for the upper limits $`N_j`$ is {2,2,2,1,1,…,1}. If we choose {$`N_j`$} to be equal to this set of minimum values, condition (4) will eliminate all but $`𝒩=3^3\times 2^8=6912`$ possible states of the system. This is still a very large number of coupled differential equations! A strong reduction in this number, however, can be obtained by applying condition (5), which we have not done previously. Of course, the number of states will be dependent on the value of $`N_{\mathrm{tot}}`$ chosen. Clearly, any meaningful value of $`N_{\mathrm{tot}}`$ must be larger than or equal to each $`N_j`$ and smaller than their sum, $`_jN_j`$. If $`N_{\mathrm{tot}}`$ is taken to be 2, the number of states will be 70 for the set of limits {2,2,2,1,1,…,1}.

To understand how this number is determined, consider the following argument. Let the number of probabilistic particles be $`m`$. For the chosen $`N_j`$ and $`N_{\mathrm{tot}}`$, there will be four different groups of states: the first group, \[0,0,…,0\], consisting of only one state {0,0,…,0}, the second group, \[1,0,…,0\], consisting of $`m`$ different states, the third group, \[1,1,0,0,…,0\], consisting of $`\frac{m(m1)}{2!}`$ states, and the fourth group, \[2,0,…,0\], consisting of only 3 states – {2,0,…,0}, {0,2,0,…,0}, and {0,0,2,0,…,0}. In our model $`m`$=11, so that the total number of states is $`1+11+55+3=70`$.

Although analytical formulae can be derived for the number of states determined by other specific values for $`m`$, {$`N_1`$,$`N_2`$,…,$`N_\mathrm{m}`$} and $`N_{\mathrm{tot}}`$, we choose to determine this number by a computer subroutine.

## 4 Results

Table 4 shows a comparison between the surface populations (in monolayers per grain) of normal (stable, non-deuterated) species calculated via the direct master equation (ME) and the Monte Carlo (MC) approaches for the low, intermediate, and high density cases. To use the former method, all $`N_j`$ were set to their minimum values except for $`N_2`$ and $`N_4`$, which correspond to O and OH. The calculations showed the latter two average abundances to be close to unity in some cases so that higher values for their upper limits had to be considered.

Table 4 also shows the values of $`N_2`$ and $`N_4`$ for the three cases as well as the limit $`N_{\mathrm{tot}}`$ for the total number of probabilistic species and the total number of states $`𝒩`$. One can see that for the high and intermediate density cases, the direct master equation method requires the solution of more than 816 coupled differential equations since the number $`𝒩=816`$ does not include the rate-like differential equations for the stable species.

For the low and intermediate density cases, the agreement between the master equation and the Monte Carlo approaches is excellent, and the master equation approach outperforms the Monte Carlo calculation in terms of computer time, dramatically so for the low density limit. For the high density case, the agreement is not as excellent, but the results of the two methods deviate by at most 10%. This small deviation can be explained by the high abundance of O and OH on the surface, a condition that requires that higher values for $`N_2`$ and $`N_4`$, and consequently many more states, be considered in the master equation approach. In general, the agreement between the two formally exact approaches can be made as good as desired by considering sufficiently high values for the $`N_j`$ and for $`N_{\mathrm{tot}}`$. This, however, will mean longer times for the master equation calculations to run, so one needs to find the right balance between the desired accuracy and the usage of computing resources. In general, it is a current weakness of the master equation approach that there is no obvious algorithm to determine the proper upper limit to the number of states to be considered.

The fractionation ratios, fXD, are defined as the ratio between the abundances of the deuterated and the normal isotopomers. The results for the fractionation ratios are given in Table 5 for the three cases discussed previously. The agreement between the two methods is reasonable but not perfect; here the major problem lies in the Monte Carlo approach, which is not accurate for small populations. Indeed, for some of the smaller fractionation ratios, no values could be calculated with this method since populations less than unity are not allowed. As an example, consider the case of CHD2OH at intermediate density. From the tables, one sees that the overall population of this molecule in the Monte Carlo approach is approximately 140 molecules. While the random (square-root) error is thus at the 10% level, the actual error is considerably greater and is probably large enough to cover the 30-40% difference between the two methods of calculation.

We have also used the direct master equation approach to model the chemistry in the temperature range 10-20 K for the high density case. Fig. 1 shows the fractionation ratio $`f`$ of doubly and multiply deuterated isotopomers in this temperature range, while a plot of the mole fraction of the normal isotopomers vs. temperature is given in Fig. 2. The parameters for the calculation (number of states, etc.) are the same as shown in Table 4. It can be seen that as the surface temperature increases from 10 to 20 K, the species (H2CO, CH3OH) made by hydrogenation involving reactions with activation energy decline sharply in abundance. This non-intuitive result derives from the fact that at higher temperatures, H atoms evaporate before they react. In addition, for most species, the fractionation ratios decline sharply with increasing temperature, especially for isotopomers that must be formed by reactions with activation energy and involving D atoms, since these atoms tunnel more poorly than their H counterparts and are more likely to evaporate.

Finally, the dependence of the fractionation ratios for doubly and multiply deuterated isotopomers on the ratio of atomic deuterium to hydrogen is shown in Fig. 3 for a 10 K cloud at high density. The abscissa is actually the ratio of the accretion rate of atomic D to that of atomic H; this is equal to the ratio of the gas-phase abundances multiplied by a factor of $`2^{0.5}`$ that derives from the relative speeds of D and H. Not surprisingly, the fractionation ratios increase as the D/H accretion ratio increases, with the slope proportional to the number of deuterium atoms in the isotopomer. In fact, the results in Fig. 3 lie very close to the very simple limit (Tielens, 1983; Brown & Millar, 1989) in which all H and D atoms that accrete onto the surface of a grain eventually react with the CO reservoir despite the activation energy barriers. In this instance, if we let $`R`$ be the ratio of the accretion rate of D to that of H, the fractionation ratios for the deuterated species are simple powers of $`R`$ multiplied by a statistical factor, where the power is just the number of D atoms in the isotopomer. The statistical factor expresses the number of possible paths leading to the isotopomer, which, for example, is 3 for CH2DOH and 1 for CH3OD. This simple limit is independent of temperature, and so works progressively more poorly as the temperature is raised from 10 K, as can be seen from the temperature dependence in Fig. 1. It also works more poorly for low and intermediate densities, where the CO reservoir does not exist.

## 5 Discussion

We have shown how a moderately complex system of surface chemical reactions can be solved by the formally exact direct master equation method. This method, or the equivalent Monte Carlo approach, is needed in the so-called accretion limit, which pertains to a situation in which the average abundance of important reactive species on a grain surface is less than unity. In the accretion limit, normal rate equations are likely to be inaccurate, and a more exact method, which treats the discrete nature of the surface populations, is needed. Of the two exact methods used to date, the direct master equation approach offers the possibility of relative ease of implementation in combined gas-grain approaches since it involves the solutions of coupled differential equations.

In order to reach the level of complexity discussed here, we have introduced a new method to limit the number of states considered. This goal is achieved by limiting the total number of reactive species considered on a grain. In the current system, which involves 11 reactive species, we have limited consideration to a total number of species in the range 3-4. Of course, this approach only works successfully if the probability for even this number of total species is small, which is what defines the accretion limit.

If one eliminates D and the five deuterium-containing radicals from the 11 reactive species considered here, and replaces them with other reactive species leading to major grain constituents, it is entirely possible that the master equation treatment described here can be extended to study combined gas-grain models of considerable complexity, in which gas-phase abundances change with time. Problems to be solved before this goal becomes a reality include the need to develop a criterion for the smallest number of states that should be considered, and the related need to determine when and for which species normal rate equations can be utilized. These problems are certain to depend on physical conditions such as grain size and temperature, as well as the rates of diffusion chosen. It would be ideal if a program could be developed to decide these issues as the integration with time proceeds.

Of the deuterated species discussed here, HDCO, D2CO, CH3OD, CH2DOH, CHD2OH, and HDO have already been detected in interstellar sources. A low-mass protostellar source – IRAS16293-2422 – contains both deuterated isotopomers of formaldehyde, both singly deuterated isotopomers of methanol, and CHD2OH in the gas phase (Parise et al., 2002). It is likely that at least some of the fractionation occurs on interstellar grains and is followed by evaporation back into the gas (Charnley et al., 1997). As shown by Parise et al. (2002); however, the granular deuteration model presented here and in Caselli et al. (2002), as well as the simpler one in Charnley et al. (1997), does not account quantitatively for all of the observations with any single value of gas-phase D chosen. Clearly, a more complex gas-grain model is needed, in which time-dependent fractionation occurs both in the gas and on dust-particle surfaces. Our current success with the direct master equation method may well allow us to consider such a complex gas-grain system in the future.

## Acknowledgments

We thank the National Science Foundation (US) for support of the Ohio State program in astrochemistry and the Ohio Supercomputer Center for time on their Cray SV1 machine. Special thanks go to V. I. Shematovich for carefully reading a version of the manuscript.


text:
# Mean-field analysis of the stability of a K-Rb Fermi-Bose mixture

## Abstract

We compare the experimental stability diagram of a Fermi-Bose mixture of 40K and 87Rb atoms with attractive interaction to the predictions of a mean-field theoretical model. We discuss how this comparison can be used to give a better estimate of the interspecies scattering length, which is currently known from collisional measurements with larger uncertainty.

The production of degenerate Fermi-Bose mixtures of atoms hulet ; salomon ; ketterle ; roati  has opened new directions in the field of ultracold atomic gases, with relevant implications for the achievement of fermionic superfluidity. Collisional interactions between bosons and fermions have a twofold role: they allow for an efficient sympathetic cooling of the fermions into quantum degeneracy, and they affect the properties of the degenerate mixture. In case of a large boson-fermion attraction, the most important predicted effects are mean-field instabilities molmer ; roth ; miya  and boson-induced interactions between fermions induced , besides a general modification of the properties of the two individual components japan ; hu ; liu . We have recently observed the collapse of the Fermi gas in a 40K-87Rb mixture science , which is indeed characterized by a large attractive interaction between the components.

In our study of the collapse we found that the order of magnitude for the critical atom numbers is in agreement with the predictions of a mean-field model for spherical geometry roth .

In this work we investigate the stability diagram of the K-Rb system, and we make a quantitative comparison of the experimental findings with mean-field theory taking into account the actual trapping potential for the atoms. We find that the mean-field model is able to reproduce the critical atom numbers for collapse using a value for the boson-fermion scattering length $`a_{BF}`$ that is in good agreement with the one we determined through the study of cold collisions simoni . Moreover, we show how the comparison of theory and experiment on the stability diagram can be used to determine the scattering length with a much lower uncertainty than that determined from collisional measurements.

The Fermi-Bose mixtures of 40K (fermions) and 87Rb (bosons) atoms are produced by sympathetic cooling in a magnetic potential roati . Both species are trapped in their low-field-seeking states with largest magnetic moment, which experience potentials of the kind

$`V_B(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{2}}m_B\left[\omega _B^2(x^2+y^2)+\omega _{Bz}^2z^2\right]`$ (1)
$`V_F(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{2}}m_F\left[\omega _F^2(x^2+(yy_0)^2)+\omega _{Fz}^2(zz_0)^2\right].`$
The trap frequencies are $`\omega _F=2\pi \times 317`$ Hz, $`\omega _{Fz}=2\pi \times 24`$ Hz for 40K, while those of 87Rb are a factor $`\sqrt{m_B/m_F}1.47`$ smaller. The system is characterized by a gravitational sag between the two clouds, which affects both the horizontal and vertical directions, due to a small misalignment of the magnetic trap with respect to the direction of gravity ribo . The values of the displacement between the two potential centers are $`y_0=3.6\mu `$m and $`z_0=20\mu `$m.

The largest samples that we are able to cool to quantum degeneracy amount to about 2$`\times `$105 bosons coexisting with 5$`\times `$104 fermions. The atom numbers for both species can be adjusted separately by controlling the loading and the evaporation procedures. The coldest Fermi gas observed have a temperature of about 0.2$`T_F`$, and the minimum temperature of the bosons that we are able to measure from the non-condensed fraction is limited to about 0.6$`T_c`$.

In the mean-field model that we use to describe this Fermi-Bose system, the ground state is determined by the solution of Gross-Pitaevskii (GP) equation for the bosons, coupled to the Thomas-Fermi equation for the fermions roth ; molmer

$`\left[{\displaystyle \frac{\mathrm{}^2}{2m_B}}^2+V_B+g_{BB}n_B+g_{BF}n_F\right]\varphi =\mu _B\varphi `$ (2)
$`n_F={\displaystyle \frac{(2m_F)}{6\pi ^2\mathrm{}^3}}^{3/2}\left(ϵ_FV_Fg_{BF}n_B\right)^{3/2}`$ (3)
where $`n_B=|\varphi ^2|`$ and $`n_F`$ are the boson and fermion densities respectively, and $`\varphi (𝐱)`$ is the condensate wave function. The boson-boson and boson-fermion interactions are described by the coupling constants $`g_{BB}=4\pi \mathrm{}^2a_{BB}/m_B`$ and $`g_{BF}=2\pi \mathrm{}^2a_{BF}/m_R`$ in terms of the $`s`$-wave scattering length $`a_{BB}`$ and $`a_{BF}`$. Here $`m_{B,F}`$ are the atomic masses and $`m_R=m_Bm_F/(m_B+m_F)`$ is the reduced mass.

In the case of the rubidium-potassium mixture, the best experimental estimates of the scattering lengths are $`a_{BB}=98.98\pm 0.02a_0`$ verhaar  and $`a_{BF}=410\pm 80a_0`$ simoni . The difference in the uncertainties for the two quantities stems from the different techniques used, which is a combination of photoassociation and Feshbach spectroscopies in the first case, and measurement of cold collisions in the second case. Although the investigation of cold collisions is the easiest way to measure the interactions in an ultracold gas, it is often accompanied by a rather large uncertainty. Indeed, what one measures in the experiment is a collisional rate, which at ultralow temperature has a dependence

$$\mathrm{\Gamma }Na^2$$
(4)
on the atom number $`N`$ and on the scattering length $`a`$. The absolute atom number, in the case of a mixture of two species, cannot be determined with an accuracy typically better that 40%, giving rise to a 20% uncertainty on the scattering length.

By solving equations (2)-(3) we can investigate the ground state and the stability of the system against collapse. The equations are solved with an imaginary-time evolution, embedded in an iterative scheme, as in roth . The stability is checked by requiring that energies, chemical potentials and central densities converge to the final value, with a precision at least of $`10^6`$. In Fig. 1 we show a typical ground state configuration, calculated with the nominal values $`a_{BF}`$=-410 $`a_0`$, $`N_B=7.510^4`$, $`N_F=210^4`$. As already discussed in Refs. molmer ; roth , the mutual attraction results in an enhancement of the density of both species in the volume of overlap. We note that the gravitational sag, that shifts the centers of the two clouds, reduces the effect of the mutual attraction with respect to the concentric case studied in Ref. roth .

The deformed ground state of the system is predicted to lead to a series of effects that can be described within a mean-field approach, including a modification of the frequencies of collective excitations japan ; hu  and of the expansion of the two clouds from the trap liu . We actually observed the latter effect already from the first experiments roati , where we had an evidence of a modified expansion of the condensate in presence of the Fermi gas.

If the atom numbers are increased above some critical values, an instability occurs. In our model, accordingly to the study reported in Refs. feldmaier ; roth  the signature of the instability is the failure of the convergence procedure during the iterative evolution toward the ground state of the system. In particular, the onset of instability is characterized by an indefinite growth of central densities which triggers the simultaneous collapse of the two species. We note that the investigation of the actual dynamics after the system has been driven into the unstable region would require a description that goes beyond Eqs. (2)-(3), in a similar fashion to what happens during the collapse of a single Bose-Einstein condensate with attractive interaction sackett ; kagan . Here we will not discuss these aspects, since we are concerned with the determination of the critical values for the onset of instabilities only.

In the experiment we observe the instability of the Fermi-Bose system as the number of condensed bosons in increased by evaporation, at a fixed number of fermions science . The signature of the instability is a sudden loss of most of the fermions from the magnetic trap, which we attribute to a largely increased three-body recombination in the collapsing samples. Here we are interested just in the critical atom numbers, but we note that further experiments, for example using a Feshbach resonance to tune $`a_{BF}`$ simoni , will be necessary to study the collapse dynamics of the system.

To compare the predictions of the mean-field model to the experimental findings on the instability, we have built a stability diagram, shown in Fig. 2. Here we report in the plane $`N_BN_F`$ the condensate and fermion atom numbers that we were able to measure in the experiment for stable samples, and compare them with the calculated critical line for instability. We note that all the data points refer to samples at temperatures of the order of 0.2-0.5$`T_F`$ for the Fermi gas, and with almost no detectable thermal fraction for the BEC. In the model, the occurrence of instability depends on both $`N_B`$, $`N_F`$, and also on the value of the scattering length $`a_{BF}`$. In Fig. 2 we have therefore plotted a family of critical lines by varying $`a_{BF}`$ around the nominal value. Note that the position of the critical line depends quite strongly on the value of $`a_{BF}`$, and it shifts toward larger atom numbers for decreasing magnitudes of $`a_{BF}`$. Although the critical line calculated for the nominal value $`a_{BF}`$=-410 $`a_0`$ is in rather good agreement with the experimental observation, the one calculated for $`a_{BF}`$=-395 $`a_0`$ better fits the data. In the experiment we observed the collapse of the Fermi gas for number pairs close to the two marked data points note , which appears to be actually the region of achievable number pairs closest to the instability.

The strong dependence on the criticality on $`a_{BF}`$ that we find in the mean-field analysis is in accordance with the scaling law for the critical number of condensate atoms

$$N_B^{crit}|a_{BF}|^\alpha $$
(5)
with $`\alpha =12`$, discussed in molmer ; miya . Here we have verified that this scaling behavior perfectly fits the numerical solutions of Eqs. (2)-(3) in case of a simplified concentric spherical configuration, by using an average trapping frequency as usually found in literature roth ; albus , as shown in Fig. 3. In this calculation we have kept the number ratio $`N_B`$/$`N_F`$ fixed, which corresponds to move along a straight line passing from the origin in the diagram of Fig. 2. This calculation also shows that the 3D geometry of the experiment results in a renormalization of the scaling exponent, $`\alpha _{3D}11`$, besides an effective increase (in modulus) of the critical scattering length, of the order $`2\%`$, due to the presence of gravity. In fact, the effect of the gravitational shift is a reduction of the density enhancements due to the attractive interaction, with respect to the case in which the to species are trapped by concentric potentials.

This analysis suggests the possibility of using a comparison of the experimental and calculated critical atom numbers to extract a precise information on the value of the scattering length for our K-Rb mixture. To account for the uncertainty in the determination of atom numbers in the experiment, in Fig. 4 we have plotted the same data points of Fig. 2, expanding or compressing both scales by 40%, and compared them with the same family of critical lines.

By taking the critical lines that best fit with the observation in these two limiting cases, we can extract the following mean-field estimate for $`a_{BF}`$

$$a_{BF}^{mf}395\pm 15a_0.$$
(6)
Note how the dependence of $`a_{BF}`$ on the number of atoms of Eq. 5, stronger than that of Eq. 4, results in a 5-fold smaller uncertainty on the scattering length determined by the stability analysis than that obtained with collisional measurements. We expect this mean-field prediction for $`a_{BF}`$ in Eq. (6) to be rather robust, since the $`40\%`$ experimental uncertainty on the number of atoms should be quite larger that the possible discrepancy between the GP theory and the actual behavior of the system, as recently discussed for the case of a single 85Rb Bose-Einstein condensate with attractive interaction bradley ; report ; savage .

The assumption of the same scaling for boson and fermion atom numbers due to experimental uncertainty is not unrealistic. Indeed, the possible sources of systematic errors in the atom number calibration are residual magnetic fields, non perfect polarization of the probe beam, or improper calibration of the imaging magnification. When using a simultaneous imaging scheme roati , all these factors affect in the same way the measured atom numbers for bosons and fermions.

We note that recently it has been proposed a beyond mean-field approach, including the second order correction in the scattering length $`a_{BF}`$, that should provide important modifications to the mean-field results albus . Nonetheless we have verified that the inclusion of this term stabilizes the system in the range of atom numbers considered here for any value of the scattering length, thus forbidding the collapse, in contrast with the experimental findings. Higher order terms might therefore be important to provide a correct beyond mean-field description, although at this stage it is not clear whether this could affect the critical number of atoms for collapse.

In conclusion, our mean-field analysis of the stability of the K-Rb Fermi-Bose mixture shows the possibility of an independent determination of the relevant scattering length of the system. The strong dependence of the critical atom numbers on the scattering length allows to largely reduce the effect of the experimental uncertainty. The value we determine for $`a_{BF}`$ is in very good accordance with that found from collisional measurements. We note that an independent, more precise determination of $`a_{BF}`$, for example using Feshbach spectroscopy simoni , would be useful to test the validity of the mean-field theory for strongly interacting Fermi-Bose mixtures, and to asses the importance of beyond mean-field effects.

###### Acknowledgements.

This work was supported by MIUR, by EC under contract HPRICT1999-00111, and by INFM, PRA “Photonmatter”. We acknowledge useful discussions with A. Albus and F. Illuminati.


text:
# 1 Introduction

## 1 Introduction

Universal amplitude ratios play a major role in modern statistical mechanics  of critical systems. Like critical indices these ratios can be used to characterize and identify universality classes. In this respect they are in general more useful than the critical indices since they are usually simpler to measure and can vary by a larger amount, thus allowing a more efficient discrimination among different universality classes.

In this letter we propose a new class of universal ratios related to the short distance expansion of the correlators of a statistical model near criticality: by construction they provide a perfect tool to probe the short distance properties of the underlying quantum field theory (QFT in the following) even in regimes where the correlation length and one point functions are not accessible in experiments or simulations.

We shall discuss the construction of these ratios with a general formalism (for a $`d`$-dimensional QFT with a generic operator content), but we will be mainly interested to two dimensional systems where the necessary assumptions are satisfied in many nontrivial cases. In particular we have in mind some integrable perturbations of Conformal Field Theories (CFT) in which case these short distance universal ratios can be evaluated exactly.

This letter is organized as follows. In sect.2, after a short reminder of known results about QFT and standard amplitude ratios, we shall construct our new universal ratios. Then in sect.3 we shall discuss, as a topical example, the case of the magnetic perturbation of the Ising model.

## 2 Universal Ratios

Let us consider a statistical model near criticality, i.e. a $`d`$dimensional QFT with operators $`\varphi _{}^{0}{}_{\alpha }{}^{}`$, regularized by a length cut-off $`a\xi `$ ($`\xi `$ being the correlation length of the system). We will assume that properties of the statistical model in the the short distance limit $`ar\xi `$ are described by a renormalized QFT obtained by an appropriate relevant perturbation $`d^dxg_i\varphi _i`$ of a CFT. (We denote generic renormalized operators by $`\varphi _\alpha `$, and their scaling dimension by $`2\mathrm{\Delta }_\alpha `$; latin indices $`i,j\mathrm{}`$ will be reserved to relevant operators, such that $`2\mathrm{\Delta }_j<d`$. As a consequence the coupling $`g_i`$ of the perturbing operator $`\varphi _i`$ will have scaling dimension $`d2\mathrm{\Delta }_i`$, such $`\xi ^1g_i^{1/(d2\mathrm{\Delta }_i)}`$.) In general any regularized operator will be equal to an appropriate mixing of renormalized operators with equal or lower scaling dimension, plus irrelevant corrections. The mixing can happen only if appropriate selection rules on scaling dimensions of operators and coupling are satisfied and symmetries of the theory are preserved. We assume that no mixing occurs among the relevant operators of the QFT to which we restrict our analysis. In this case we are left with simple multiplicative relations (up to irrelevant corrections): $`\varphi _{}^{0}{}_{j}{}^{}=Z_j^1\varphi _j+\mathrm{}`$ in which the factors $`Z_j`$ are the only non-universal avatar of the original statistical model. The invariance of the perturbation, $`g_i^0\varphi _{}^{0}{}_{i}{}^{}=g_i\varphi _i+\mathrm{}`$, implies the relation $`g_i^0=Z_ig_i`$. It is clear that any scale invariant constant built up from relevant operators and coupling that is invariant by the rescaling

$`\varphi _jK_j^1\varphi _j;g_iK_ig_i`$ (1)
will not depend on $`Z_j`$ and as a consequence will be universal. (We borrow from  the name “metric factors” for our normalization constants $`K_j`$.)

Let us now review some known universal amplitude ratios related to VEV and correlators (For a thorough review and additional references see .)

Let us start with ratios involving one point functions of relevant operators and denote

$$\varphi _j_i=A_j^ig_i^{\frac{2\mathrm{\Delta }_j}{d2\mathrm{\Delta }_i}}$$
(2)
the vacuum expectation value (VEV) of the operator $`\varphi _j`$ in the theory obtained perturbing with respect to $`\varphi _i`$. Clearly, the scaling Eq.(1) is equivalent to replace

$$A_j^iK_jA_j^i(K_i)^{\frac{2\mathrm{\Delta }_j}{d2\mathrm{\Delta }_i}}$$
(3)
in Eq.(2). The only remaining task to obtain universal ratios is to construct suitable combinations of the VEV (or of their derivatives wrt $`g_i`$) so as to eliminate all the metric factors (see  for recent applications). An exhaustive list of these combinations can be found, for instance, in .

A similar approach can be followed also if we take, instead of one point functions, the expectation value of two (or more points) functions. Much work has been done in the past years in the context of $`O(N)`$ scalar models: for instance it was soon realized (see  and references therein) that the connected spin spin correlator as a function of $`r/\xi `$ is a universal function (up to an overall constant) whose asymptotics can be parametrized in terms of universal constants.

When studying connected correlators in the regime $`r\xi a`$ one has access to universal ratios of the so called “overlap amplitudes”. This line of research have been recently pushed forward in the case of the 2d Ising model perturbed by a magnetic field  and for the thermal perturbation of the 3d Ising model .

The complementary possibility is to look at the short distance behavior $`ar\xi `$ of correlators. In the context of $`O(N)`$ scalar models Wilson’s Operator Product Expansion (OPE) revealed to be  the appropriate theoretical tool to deal with short distance regime and $`ϵ=4d`$ expansion techniques have been applied to describe correlators’ asymptotics and to compare them with simulations or with experiments (see ).

In this short letter we propose a new class of universal amplitudes –obtained from ratios of correlators in short distance regime– that do not require additional measures of $`\xi `$ or VEV: the idea is that if one is looking at very small perturbations of the critical point it might happen that the correlation length and VEV estimates are not accessible in the particular experimental or numerical realization in which one is interested.

We pursue to describe short distance asymptotics by applying a specialization of the OPE technique to perturbed CFT recently proposed in  and known as InfraRed Safe approach (IRS). Roughly speaking, the main ingredient of the IRS approach is simply to choose a renormalization scheme in which the Wilson’s coefficients of the OPE have a regular, IR safe, perturbative expansion with respect to the coupling and then to expand them in OPE. (See  for all-orders IRS formulas and many cautionary remarks on their range of validity.) The approach can be extended in principle to $`n`$points correlators and irrelevant operators, but in the following we shall concentrate on the two point functions of relevant operators only, because the hypotheses of non mixing and of absence of divergences –see below– are likely to be satisfied in this case and because in more complicate generalizations the amplitude ratios will be too difficult to study from a numerical point of view.

We choose as expansion variable the scaling combination $`tg_ir^{d2\mathrm{\Delta }_i}(r/\xi )^{d2\mathrm{\Delta }_i}`$ and rescale the correlators as follows

$`F_{jk}^i(t)r^{2\mathrm{\Delta }_j+2\mathrm{\Delta }_k}\varphi _j\varphi _k_i.`$ (4)
(Notice that $`\varphi _j\varphi _k_i`$ is just the expectation value of the two operators and  not the connected correlator.) In this notations IRS gives

$`F_{jk}^i(t){\displaystyle \underset{\alpha ,n}{}}A_\alpha ^i\widehat{{\displaystyle \frac{^nC_{jk}^\alpha }{g_i^n}}}{\displaystyle \frac{t^{\frac{2\mathrm{\Delta }_\alpha }{d2\mathrm{\Delta }_i}+n}}{n!}}{\displaystyle \underset{a}{}}f_{jka}^it^{z_{jka}^i}`$ (5)
where the constants $`f_{jka}^i`$ are suitable combinations of $`\widehat{^nC_{jk}^\alpha /g_i^n}`$ (derivatives of Wilson’s coefficients evaluated at $`g_i=0`$ and scaled by appropriate powers of $`r`$ to obtain adimensional quantities) and of the amplitudes $`A_\alpha ^i`$ of operators’ VEV, defined as in Eq.(2); $`z_{jka}^i`$ are appropriate combinations of scaling dimensions. Notice that in Eq.(5) we assumed a powerlike short-distance expansion in terms of the variable $`t`$ (at least up to the order we will need): in general it may happen that the underlying QFT contains pairs of “resonant” operators and terms proportional to powers of $`\mathrm{log}t`$ may appear in the correlators at some order in the expansion (the typical example of this situation is the thermal perturbation of the 2d Ising model). If logs appear in the terms of the expansion the present treatment cannot be applied as it is (we shall deal with these peculiar cases in forthcoming publications).

The key point for the construction of universal ratios is that the scaling Eq.(1) is equivalent to the nontrivial rescaling

$$f_{jka}^iK_jK_kK_i^{z_{jka}^i}f_{jka}^i$$
(6)
(mixing factors from $`\varphi _\alpha `$ are compensated by those which come from the corresponding redefinition of Wilson’s coefficients, even in the general case of matricial mixing). All we need at this point is to look for combinations in which the non universal factors cancel among each other as in the case of standard amplitude ratios. This can be performed in two steps. First by combining the correlators in appropriate ratios to eliminate the trivial overall dependence on $`K_j,K_k`$ 111Depending on the particular case under study, it may be simpler to eliminate this dependence by using one point functions. However in this letter, following the discussion at the beginning of this section, we shall assume  not to have access to these quantities., and then by choosing a combination of ratios to have a finite non-zero result in the $`t0`$ limit (this step replaces the knowledge of the correlation length in the standard approach).

There are in general several ways to construct these combinations, but it is clear that the most interesting ones are those which involve the smallest number of different operators and thus have greater chances to be observed in experiments or in simulations. In this letter we propose two options which –in our opinion– are the best candidates. Let us look at them in detail.

### 2.1 QFT’s with three or more relevant operators.

Let us choose a model (like for instance the tricritical Ising model) which contains at least three non trivial relevant operators (it doesn’t matter if one of the three is the same $`\varphi _i`$ which we use as a perturbation). Let us call them $`\varphi _1`$,$`\varphi _2`$ and $`\varphi _3`$, and let us construct the following combinations:

$$R_{12}^i\frac{\varphi _1\varphi _1_i\varphi _2\varphi _2_i}{[\varphi _1\varphi _2_i]^2}=\frac{F_{11}^iF_{22}^i}{[F_{12}^i]^2}$$
(7)
and similarly for the pairs ($`\varphi _1,\varphi _3`$) and ($`\varphi _2,\varphi _3`$). The small $`t`$ expansion of the $`R_{jk}^i`$ can be reconstructed from that of their components and one finds in general

$$R_{jk}^ir_{jk}^it^{z_{jk}^i}+\mathrm{}$$
(8)
where $`r_{jk}^i`$ and $`z_{jk}^i`$ are simple combinations of the $`f_{jka}^i`$ and $`z_{jka}^i`$ which appear in the small $`t`$ expansions (see Eq.(5)) of the correlators involved in the ratio.

The reader can easily check that the overall metric factors $`K_j,K_k`$ disappear in each of the combinations $`R_{jk}^i`$. At this point one can construct one or more universal combinations out of these ratios by combining the $`R_{jk}^i`$ so as to cancel the leading $`t`$ dependence. The $`t0`$ limit of any one of these combinations is clearly an universal quantity and can be exactly expressed in terms of the Wilson’s coefficients, their derivatives and the one point functions of the theory. For instance one such combination could be:

$$Q_{123}^i\underset{t0}{lim}\frac{(R_{12}^i)^{z_{13}^i}}{(R_{13}^i)^{z_{12}^i}}=\frac{(r_{12}^i)^{z_{13}^i}}{(r_{13}^i)^{z_{12}^i}}$$
(9)
### 2.2 Subtracted correlators.

Unfortunately the above construction cannot be applied, for instance, to the very interesting case of the Ising model which has only two non trivial relevant operators in the spectrum. However there is a simple modification which allows to overcome this problem.

Let us construct the differences:

$$F_{jj,s}^i(t)F_{jj}^i(t)F_{jj}^i(0)$$
(10)
which we shall call in the following “subtracted correlators”. Since $`F_{jj}^i(0)`$ is a constant related to the correlator at the critical point, the subtracted correlators are almost as easily accessible from a numerical point of view as the ordinary correlators. With the help of these subtracted correlators we can construct new combinations, like for instance

$$R_{j,s}^i\frac{\varphi _j\varphi _j_{i,s}}{\varphi _j\varphi _j_i}(j=1,2)$$
(11)
which can be suitably combined among them or with $`R_{12}^i`$ so as to eliminate the residual $`t`$ dependence in the $`t0`$ limit, thus leading to a universal combination.

### 2.3 QFT estimates.

Let us finally address the question of the theoretical estimate of these new universal amplitude ratios. In the case of three dimensional $`O(N)`$ scalar models it should be possible to obtain $`ϵ`$-expansion (or fixed dimensional) estimates by following the same lines used to describe short distance behavior of field field correlator (see ). The situation is much better in two dimensions where, thanks to the progress of CFT one knows exactly the operator content, scaling dimensions and correlators of the underlying critical QFT. Moreover IRS method gives an exact all order integral representation for the derivatives of Wilson’s coefficients $`\widehat{^nC_{jk}^\alpha /g_i^n}`$ in terms of CFT correlators from which accurate or even exact predictions can be obtained. For what concerns VEV amplitudes, it is remarkable that the approximate technique known as Truncated Conformal Space allows for theoretical predictions of $`A_\alpha ^i`$ (see  and references therein) and in addition in the case of 2d integrable perturbations many $`A_\alpha ^i`$ can be evaluated exactly . Notice that in 2d, integrable perturbations exist that are not described in terms of $`O(N)`$ models.

## 3 Example: the magnetic perturbation of the Ising model

The Ising model perturbed by the spin operator $`\sigma `$ (in $`d=2,3`$) is the simplest non-trivial example of the previous discussion. The perturbing parameter in this case is the magnetic field $`h`$ (to simplify notations we will drop in this section the perturbation label from amplitudes and ratios). The vacuum expectation values of the two relevant operators are:

$`\sigma _h=A_\sigma h^{\frac{2\mathrm{\Delta }_\sigma }{d2\mathrm{\Delta }_\sigma }};ϵ_h=A_ϵh^{\frac{2\mathrm{\Delta }_ϵ}{d2\mathrm{\Delta }_\sigma }}`$ (12)
In the following we will concentrate on the construction of our new short distance amplitude ratios; for different amplitude ratios testing other regimes of the model, see . The first few terms of the short distance expansions for the correlators look like:

$`F_{\sigma \sigma }(t)`$ $``$ $`\widehat{C_{\sigma \sigma }^\mathrm{𝟏}}+A_ϵ\widehat{C_{\sigma \sigma }^ϵ}t^{\frac{2\mathrm{\Delta }_ϵ}{d2\mathrm{\Delta }_\sigma }}+A_\sigma \widehat{_hC_{\sigma \sigma }^\sigma }t^{1+\frac{2\mathrm{\Delta }_\sigma }{d2\mathrm{\Delta }_\sigma }}`$
$`F_{ϵϵ}(t)`$ $``$ $`\widehat{C_{ϵϵ}^\mathrm{𝟏}}+A_ϵ\widehat{C_{ϵϵ}^ϵ}t^{\frac{2\mathrm{\Delta }_ϵ}{d2\mathrm{\Delta }_\sigma }}+A_\sigma \widehat{_hC_{ϵϵ}^\sigma }t^{1+\frac{2\mathrm{\Delta }_\sigma }{d2\mathrm{\Delta }_\sigma }}`$
$`F_{\sigma ϵ}(t)`$ $``$ $`A_\sigma \widehat{C_{\sigma ϵ}^\sigma }t^{\frac{2\mathrm{\Delta }_\sigma }{d2\mathrm{\Delta }_\sigma }}+\widehat{_hC_{\sigma ϵ}^\mathrm{𝟏}}t`$
where $`thr^{d2\mathrm{\Delta }_\sigma }`$, $`\widehat{_h^nC_{ij}^k}`$ are the rescaled derivatives w.r.t. $`h`$ of Wilson’s coefficients. We warn the reader that the ordering of powers of $`t`$ from analogous contributions is in general different in $`d=3`$ and $`d=2`$ due to different scale dimensions. The exact values of the first terms of short distance expansion for the case $`d=2`$ are discussed in ; here we only need to know that in $`d=2`$ the coefficients $`\widehat{C_{ϵϵ}^ϵ}`$ and $`\widehat{_hC_{ϵϵ}^\sigma }`$ are exactly zero and $`F_{ϵϵ,s,d=2}=\frac{1}{2}\widehat{_h^2C_{ϵϵ}^\mathrm{𝟏}}t^2+O(t^{32/15})`$ (see ). Let us now introduce the correlator ratios as in Eq.(7) and Eq.(11). (Note that leading behavior of subtracted operators turns out to be dimension dependent.)

$`R_{\sigma ϵ}(t){\displaystyle \frac{F_{\sigma \sigma }(t)F_{ϵϵ}(t)}{F_{\sigma ϵ}^2(t)}}{\displaystyle \frac{\widehat{C_{\sigma \sigma }^\mathrm{𝟏}}\widehat{C_{ϵϵ}^\mathrm{𝟏}}}{(\widehat{C_{\sigma ϵ}^\sigma }A_\sigma )^2}}t^{\frac{4\mathrm{\Delta }_\sigma }{d2\mathrm{\Delta }_\sigma }}+\mathrm{}`$
$`R_{\sigma ,s}(t){\displaystyle \frac{F_{\sigma \sigma ,s}(t)}{F_{\sigma \sigma }(t)}}A_ϵ{\displaystyle \frac{\widehat{C_{\sigma \sigma }^ϵ}}{\widehat{C_{\sigma \sigma }^\mathrm{𝟏}}}}t^{\frac{2\mathrm{\Delta }_ϵ}{d2\mathrm{\Delta }_\sigma }}+\mathrm{}`$
$`R_{ϵ,s,d=2}(t){\displaystyle \frac{F_{ϵϵ,s}(t)}{F_{ϵϵ}(t)}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{_h^2\widehat{C_{ϵϵ}^\mathrm{𝟏}}}{\widehat{C_{ϵϵ}^\mathrm{𝟏}}}}t^2+\mathrm{}`$
$`R_{ϵ,s,d=3}(t){\displaystyle \frac{F_{ϵϵ,s}(t)}{F_{ϵϵ}(t)}}A_ϵ{\displaystyle \frac{\widehat{C_{ϵϵ}^ϵ}}{\widehat{C_{ϵϵ}^\mathrm{𝟏}}}}t^{\frac{2\mathrm{\Delta }_ϵ}{32\mathrm{\Delta }_\sigma }}+\mathrm{}`$
It is easy to see that the simplest combinations of these ratios which in the $`t0`$ limit lead to universal quantities are

$`Q_{\sigma ϵ,s}(t)`$ $``$ $`R_{\sigma ϵ}^{\frac{\mathrm{\Delta }_ϵ}{2\mathrm{\Delta }_\sigma }}(t)R_{\sigma ,s}(t)`$ (13)
$`Q_{\sigma ϵ,s,d=2}^{^{}}(t)`$ $``$ $`R_{ϵ,s,d=2}^{\frac{4}{15}}(t)R_{\sigma ,s}(t)`$ (14)
$`Q_{\sigma ϵ,s,d=3}^{^{}}(t)`$ $``$ $`R_{ϵ,s,d=3}^1(t)R_{\sigma ,s}(t)`$ (15)
$`Q_{\sigma ϵ,s}^{^{}}`$ is probably the best candidate in $`d=3`$, while in $`d=2`$ the best choice is to study $`Q_{\sigma ϵ,s}`$. Let us concentrate on $`d=2`$ case, where $`\frac{\mathrm{\Delta }_ϵ}{2\mathrm{\Delta }_\sigma }=4`$, and the universal ratio can be written directly in terms of correlators as:

$$Q_{\sigma ϵ,s,d=2}(t)=\frac{\sigma \sigma _{h,s}(\sigma \sigma _h)^3(ϵϵ_h)^4}{(\sigma ϵ_h)^8}.$$
(16)
The ratio $`Q_{\sigma ϵ,s,d=2}(0)`$ can be obtained exactly. If we choose the CFT field normalizations such $`\widehat{C_{\sigma \sigma }^\mathrm{𝟏}}=\widehat{C_{ϵϵ}^\mathrm{𝟏}}=1`$ we have $`\widehat{C_{\sigma ϵ}^\sigma }=\widehat{C_{\sigma \sigma }^ϵ}=\frac{1}{2}`$. Thus we find:

$$Q_{\sigma ϵ,s,d=2}(0)=128\frac{A_ϵ}{A_{\sigma }^{}{}_{}{}^{8}}$$
(17)
By using the exact values of the VEV amplitudes (see  and references therein)

$$A_\sigma =1.27758227\mathrm{}A_ϵ=2.00314\mathrm{}$$
(18)
we finally obtain $`Q_{\sigma ϵ,s,d=2}(0)=36.125\mathrm{}`$.

Acknowledgements: We thank A. Pelissetto and J. Zinn–Justin for helpful suggestions. This work was partially supported by the European Commission TMR program HPRN-CT-2002-00325 (EUCLID). The work of one of us (P.G.) is supported by the COFIN ”Teoria dei campi, meccanica statistica e sistemi elettronici”.