text:
# Hydrodynamical simulations of the Sunyaev–Zel’dovich effect: cluster scaling relations and X-ray properties

## 1 Introduction

The prospect of surveying for galaxy clusters across a wide range of redshifts using the Sunyaev–Zel’dovich (SZ) effect (Sunyaev & Zel’dovich 1972, 1980) raises the hope of a much-improved understanding of cluster physics, and in particular of the evolution of the cluster gas. The SZ effect is a particularly powerful tool when combined with observations of X-ray emission from the hot cluster gas, as the two techniques probe different properties of the cluster gas distribution, and with the XMM-Newton and Chandra X-ray satellites in their main operations phases the number of clusters and groups with quality X-ray observations is rapidly increasing. It is therefore timely to extend theoretical models of clusters, where efforts have hitherto been concentrated upon their X-ray properties, to incorporate predictions for the SZ effect.

Numerical simulations indicate that the distribution of non-baryonic dark matter in clusters, which dominates their mass, is approximately self-similar (e.g. Navarro, Frenk & White 1995, 1997). Observations indicate, however, that the baryonic gas component cannot share such a degree of self-similarity. This is particularly evident from observations of the $`L_\mathrm{X}`$$`T_\mathrm{X}`$ relation, which show it to be steeper than predicted by the self-similar model (e.g. Edge & Stewart 1991; Xue & Wu 2000). The reason for this discrepancy is that the gas is less centrally concentrated than the dark matter, due to physical processes (in addition to gravity) that raised the entropy of the gas (e.g. Evrard & Henry 1991; Kaiser 1991; Bower 1997; Voit et al. 2002; Ponman, Sanderson & Finoguenov 2003). This entropy could come from direct heating of the hot gas from stars or AGN (e.g. Wu, Fabian & Nulsen 2000; Bower et al. 2001; Quilis, Bower & Balogh 2001), or from radiative cooling which transforms low-entropy material into stars allowing high-entropy material to flow in to replace it (e.g. Knight & Ponman 1997; Pearce et al. 2000; Bryan 2000). Understanding the relative importance of heating and cooling, and how these processes conspire to establish excess entropy in clusters, is one of the key problems in cluster gas physics.

The implications of an excess of entropy on the global properties of the SZ effect (e.g. source counts, the mean Compton parameter and the angular power spectrum) has previously been investigated using both analytical techniques (e.g. Cavaliere & Menci 2001; Holder & Carlstrom 2001; Zhang & Wu 2003) and numerical simulations (e.g. Springel, White & Hernquist 2001; da Silva et al. 2001; White, Hernquist & Springel 2002). The mass dependence of the SZ effect in clusters has been addressed by Metzler (1998) and by White et al. (2002). More recently, McCarthy et al. (2003a) used a semi-analytic model to derive scaling relations between the central Compton parameter and mass, temperature, SZ flux density and X-ray luminosity. These scalings were subsequently used in a companion paper (McCarthy et al. 2003b) to estimate the level of the entropy floor compatible with present-day SZ observations.

The purpose of this paper is to investigate the relationship between SZ and X-ray properties of clusters in $`N`$-body/hydrodynamical simulations that include models for both cooling and preheating. The simulations we use have already been demonstrated to give good overall agreement with observed X-ray scaling relations at redshift zero (Muanwong et al. 2001; Thomas et al. 2002; Muanwong et al. 2002). Our approach is mainly focused on the derivation of theoretical scaling laws, and their evolution with redshift, between the integrated SZ flux density and other cluster properties, such as X-ray temperature and luminosity. We quantify deviations from self-similar evolution on both the slope and normalization of the scalings and provide fits which can be used for comparison with other theoretical models and observations. An observationally-motivated analysis of SZ/X-ray correlations, that can be more readily applied to ‘blind’ searches for high-redshift clusters, will be presented in a future paper.

This paper is organized as follows. We define the observational quantities used and the simple scaling laws predicted by the self-similar model in Section 2. Details of our simulations and how our cluster catalogues were constructed are presented in Section 3. In Section 4, we investigate the correlation between thermal SZ integrated flux and other 3D cluster properties (mass, mass-weighted and emission-weighted temperature and X-ray luminosity) from the simulation data at redshift zero. The evolution of these relations with redshift is studied in Section 5, before we draw conclusions in Section 6.

## 2 Theoretical Framework

### 2.1 Definitions of physical quantities

The SZ effect arises due to inverse Compton scattering of CMB photons off free electrons, with the largest SZ signal being produced by intracluster gas. The total SZ flux density produced by a cluster is the integral of the SZ sky brightness over its subtended solid angle

$$S_\nu =I_0\left[g(x)yh(x)b\right]𝑑\mathrm{\Omega },$$
(1)
where $`I_0=2\left(k_\mathrm{B}T_0\right)^3/(hc)^22.28\times 10^4`$ mJy arcmin-2, 1111 mJy = $`10^{26}`$ erg s-1 m-2 Hz-1. for a mean CMB temperature of $`T_02.725`$ K (Mather et al. 1999). The first term inside the brackets accounts for the thermal SZ effect, due to internal motion of the electrons, whereas the second term, the kinetic SZ effect, is the contribution due to bulk motion of the gas. We ignore the kinetic SZ effect in this paper as the thermal SZ effect dominates at all frequencies except $`\nu 217`$GHz ($`x=h\nu /k_\mathrm{B}T_03.83`$).

The Comptonization parameter $`y`$ contains information on the structure of the intracluster gas

$$y=\frac{k_\mathrm{B}\sigma _\mathrm{T}}{m_\mathrm{e}c^2}T_\mathrm{e}n_\mathrm{e}𝑑l,$$
(2)
where $`T_\mathrm{e}`$ and $`n_\mathrm{e}`$ are the temperature and density of the electrons, $`\sigma _\mathrm{T}6.65\times 10^{25}\mathrm{cm}^2`$ the Thomson cross-section, $`c`$ the speed of light and $`m_\mathrm{e}`$ the electron rest mass. It is therefore convenient to define the integrated contribution to $`S_\nu `$ from $`y`$

$$Y=y𝑑\mathrm{\Omega }=\frac{y}{d_A^2}𝑑A=\frac{k_B\sigma _\mathrm{T}}{m_\mathrm{e}c^2d_A^2}_VT_\mathrm{e}n_\mathrm{e}𝑑V,$$
(3)
where $`d_A`$ is the angular diameter distance and $`dA=d\mathrm{\Omega }d_A^2`$ is the sky projected area element of the cluster. Since $`Y`$ depends on distance, we will usually work with the intrinsic thermal SZ signal, defined as $`Y^{\mathrm{int}}=Yd_A^2(z)`$.

The X-ray luminosity of a cluster is

$$L_\mathrm{X}=n_\mathrm{e}n_\mathrm{H}\mathrm{\Lambda }_\mathrm{c}(T)𝑑V.$$
(4)
The integral is over the volume of the cluster and $`\mathrm{\Lambda }_\mathrm{c}(T)`$ denotes the cooling rate of the gas, assuming collisional ionization equilibrium. At temperatures above a few keV, the emission is dominated by thermal bremsstrahlung, but below this line emission becomes important.

### 2.2 Self-similar scaling relations

In the absence of any non-gravitational heating and cooling processes, clusters scale self-similarly to a good approximation (Kaiser 1986; Navarro et al. 1995). In self-similar models, the temperature of the gas scales with the cluster mass as

$$TM^{2/3}(1+z),$$
(5)
assuming the density of the gas scales with the mean density $`\rho (1+z)^3`$ and the system is in virial equilibrium. Applying this to Eq. (3) we get

$$Y^{\mathrm{int}}\{\begin{array}{cc}f_{\mathrm{gas}}T^{5/2}(1+z)^{3/2}\hfill & \\ f_{\mathrm{gas}}M^{5/3}(1+z),\hfill & \end{array}$$
(6)
where $`f_{\mathrm{gas}}`$ is the gas mass fraction of the cluster.

The X-ray luminosity of a self-similar cluster scales as

$$L_\mathrm{X}\{\begin{array}{cc}f_{\mathrm{gas}}^2T^2(1+z)^{3/2}\hfill & \\ f_{\mathrm{gas}}^2M^{4/3}(1+z)^{7/2},\hfill & \end{array}$$
(7)
assuming that the X-rays are due to bremsstrahlung emission ($`\mathrm{\Lambda }_\mathrm{c}T^{1/2}`$).

## 3 Method

### 3.1 Simulation details

Results are presented from three simulations, described in detail elsewhere (Thomas et al. 2002; Muanwong et al. 2002); we summarize pertinent details here.

We assume a $`\mathrm{\Lambda }`$CDM cosmology, setting the density parameter, $`\mathrm{\Omega }_\mathrm{m}=0.35`$, cosmological constant, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$, hubble parameter, $`h=0.71`$, baryon density, $`\mathrm{\Omega }_bh^2=0.019`$, CDM power spectrum shape parameter, $`\mathrm{\Gamma }=0.21`$, and normalization $`\sigma _8=0.9`$. These values are similar to those favoured by current observations including WMAP. Initial conditions were generated using $`160^3`$ particles each of baryonic and dark matter, perturbed from a regular grid of comoving length, $`L=100h^1\mathrm{Mpc}`$. These choices set the gas and dark matter particle masses equal to $`2.6\times 10^9h^1\mathrm{M}_{}`$ and $`2.1\times 10^{10}h^1\mathrm{M}_{}`$ respectively.

The runs were started at redshift $`z=49`$ and evolved to $`z=0`$ using a parallel version of the hydra $`N`$-body/hydrodynamics code (Couchman, Thomas & Pearce 1995; Pearce & Couchman 1997), which uses a Smoothed Particle Hydrodynamical (SPH) implementation similar to that discussed in Thacker & Couchman (2000), and conserves both energy and entropy to a satisfactory degree. The gravitational softening was fixed at $`50h^1\mathrm{kpc}`$ in comoving co-ordinates until $`z=1`$, then at $`25h^1\mathrm{kpc}`$ in physical co-ordinates until $`z=0`$.

The simulations differed only in the adopted heating/cooling model. In the first model, a Non-radiative simulation, the gas was subject only to gravitational, adiabatic and viscous forces, and so could only change its entropy through shock-heating. This model fails to reproduce the X-ray scaling relations (Muanwong et al. 2001) but nevertheless can be used to test the self-similar scalings presented in the previous section.

The second simulation, a Cooling model, allowed gas to cool radiatively as described by Thomas & Couchman (1992). Tabulated values of $`\mathrm{\Lambda }_\mathrm{c}`$ were generated from cooling tables in Sutherland & Dopita (1993), assuming a time-varying global gas metallicity, $`Z=0.3(t/t_0)Z_{}`$, where $`t_012.8\mathrm{Gyr}`$ is the current age of the Universe. Cooled material (identified as all gas particles with temperatures $`T<1.2\times 10^4`$ K and overdensities $`\delta >1000`$) was converted into collisionless particles. Note that the choice of density threshold (considerably lower than the observed density of star formation) largely circumvents the problem of unphysical cooling of hot gas particles in contact with cold dense gas, as well as increasing run-time efficiency. The global fraction of cooled gas in the simulation box is $``$15 per cent at $`z=0`$. The removal of this gas from the hot phase sufficiently increases the entropy in the centres of clusters to reproduce the X-ray scaling relations, but the cooled fraction is 2 to 3 times higher than observed (Cole et al. 2001). In this respect, this run provides a maximum variation of $`Y`$ from a cluster due to effect of radiative cooling alone.

Finally, the third simulation, a Preheating model, also allowed the gas to cool radiatively, but the gas was impulsively heated by 1.5 keV per particle at $`z=4`$. This generates the required minimum entropy in the centres of clusters but prevents further cooling to $`10^4`$K; the global fraction of cooled gas in the simulation box is only $``$0.5 per cent at $`z=0`$.

Note that we do not attempt to perform a numerical resolution study for this work. Apart from being computationally demanding, it would be misleading to study the effects of resolution on our current models. In particular, increasing the resolution in the Cooling model would result in an even higher cooled fraction, and so attaining convergence in this model is undesirable. Instead we interpret our models as a relatively crude first attempt at capturing the effects of increasing the entropy of the gas by the required amount in two completely different ways (i.e. from cooling and heating). A proper convergence analysis is warranted when feedback is included (e.g. Springel & Hernquist 2003) and will be investigated in future work.

### 3.2 Cluster selection and SZ/X-ray estimators

We identified clusters in our simulations according to the method described in Thomas et al (1998). Clusters were selected by first creating a minimal-spanning tree of all dark matter particles whose density exceeds $`\delta =178\mathrm{\Omega }^{0.55}(z)`$ times the mean dark matter density (the approximate value for a virialized sphere, as predicted by the spherical top-hat model; Eke, Navarro & Frenk 1998). The tree was then pruned into clumps using a maximum linking length equal to $`0.5\delta ^{1/3}`$ times the mean interparticle separation. A sphere was then grown around the densest particle in each clump until the enclosed mean density exceeded a value $`\mathrm{\Delta }`$ in units of the comoving critical density. For this paper, we use $`\mathrm{\Delta }=200`$, larger than the virial value for our cosmology ($`\mathrm{\Delta }111`$) but commonly used by other authors. Master cluster catalogues222The catalogues are available to download at the website virgo.sussex.ac.uk were produced containing only objects with at least 500 particles each of gas and dark matter, equivalent to a mass limit, $`M_{\mathrm{lim}}1.18\times 10^{13}h^1M\mathrm{}`$. At $`z=0`$, the catalogues contain 428, 457 and 405 clusters for the Non-radiative, Cooling and Preheating simulations respectively.

For each cluster, we then estimated various observable quantities: the intrinsic SZ signal, $`Y^{\mathrm{int}}`$, the X-ray emission-weighted gas temperature, $`T_\mathrm{X}`$ and the X-ray luminosity, $`L_\mathrm{X}`$. We identified the hot intracluster gas as all gas particles within $`R_{200}`$ with $`T>10^5`$K (in practice there are very few particles with $`10^4`$K$`<T<10^5`$K as their cooling times are very short) and assume full ionization, such that the number of electrons per baryon is $`\eta =0.88`$ for a hydrogen mass fraction $`X=0.76`$. We have

$$Y^{\mathrm{int}}=\frac{k_\mathrm{B}\sigma _\mathrm{T}}{m_\mathrm{e}c^2}\frac{\eta }{m_\mathrm{p}}f_{\mathrm{gas}}MT_{\mathrm{mw}},$$
(8)
where $`T_{\mathrm{mw}}=m_iT_i/m_i`$ is the mass-weighted temperature of the gas, $`f_{\mathrm{gas}}=m_i/M`$ the gas fraction, $`m_i`$ the mass of hot gas particle $`i`$ and $`M`$ the total mass of the cluster.

The X-ray emission-weighted temperature was estimated as

$$T_\mathrm{X}=\frac{\mathrm{\Sigma }_im_i\rho _i\mathrm{\Lambda }_{\mathrm{soft}}(T_i,Z)T_i}{\mathrm{\Sigma }_im_i\rho _i\mathrm{\Lambda }_{\mathrm{soft}}(T_i,Z)},$$
(9)
where $`\rho _i`$ and $`T_i`$ are the density and temperature of the hot gas particles. Since most observed temperatures use instruments sensitive to soft X-rays, we use a soft-band cooling function, $`\mathrm{\Lambda }_{\mathrm{soft}}`$, from Raymond & Smith (1977) for an energy range 0.3–1.5 keV.

The bolometric X-ray luminosity of each cluster, corrected from its soft-band emission, was estimated as

$$L_\mathrm{X}=\frac{\mathrm{\Lambda }_{\mathrm{bol}}(T_\mathrm{X})}{\mathrm{\Lambda }_{\mathrm{soft}}(T_\mathrm{X})}\underset{i}{}\frac{m_i\rho _i\mathrm{\Lambda }_{\mathrm{soft}}(T_i,Z)}{(\mu m_\mathrm{H})^2},$$
(10)
where $`\mu m_\mathrm{H}=10^{24}`$g is the mean molecular mass of the gas and $`\mathrm{\Lambda }_{\mathrm{bol}}`$ is the bolometric cooling function (Sutherland & Dopita 1993).

Finally, to avoid introducing artificial trends in the derived scalings with temperature and X-ray luminosity (Sections 4 & 5), we imposed lower limits on $`T_\mathrm{X}`$ and $`L_\mathrm{X}`$ (by inspecting the $`T_\mathrm{X}`$$`M`$ and $`L_\mathrm{X}`$$`M`$ relations at all redshifts) to create (separate) catalogues that were complete in the two quantities. At each redshift, these limits approximately correspond to the maximum temperature and luminosity of all clusters at our mass limit, $`M_{\mathrm{lim}}`$. For example, when investigating temperature scalings at $`z=0`$, all clusters with $`T_\mathrm{X}<0.35,0.65,0.60`$ keV (for the Non-radiative, Cooling and Preheating simulations respectively) were discarded, reducing the number of clusters in our catalogues by around 30 per cent in the Cooling and Preheating simulations and 46 per cent in the Non-radiative run. Similarly, for scalings with X-ray luminosity at $`z=0`$, our original cluster catalogues were trimmed by selecting only clusters with $`L_\mathrm{X}>2.9\times 10^{43}`$, $`9.0\times 10^{41}`$, $`4.5\times 10^{41}\mathrm{erg}\mathrm{s}^1`$, for the Non-radiative, Cooling and Preheating simulations, reducing the number of clusters in the original catalogues by about 52, 46, and 29 per cent respectively.

## 4 Scaling relations at redshift zero

### 4.1 The $`Y`$$`M`$ and $`Y`$$`T_{\mathrm{mw}}`$ relations

We begin by correlating $`Y^{\mathrm{int}}`$ with mass and mass-weighted temperature, all measured within $`R_{200}`$. An early study of the relation between mass and SZ absorption in clusters was made in unpublished work by Metzler (1998), who verified with non-radiative simulations the related self-similar scaling, $`yM`$. More recently, White et al. (2002) also recovered the self-similar scaling, $`YM^{5/3}`$, from simulations with the same box size but a slightly higher (factor of 2.5) mass resolution than that used here, concluding that the combined effect of cooling and feedback (galactic winds) on the SZ properties of their clusters was small.

Fig. 1 illustrates our $`Y_{200}^{\mathrm{int}}M_{200}`$ and $`YT_{\mathrm{mw}}`$ relations for clusters in the Non-radiative (crosses) and Cooling (circles) simulations. The solid lines represent power-law best fits to the data, while the dashed line is the best-fit to clusters in the Preheating simulation, omitted for the purpose of clarity. The figure shows that $`Y_{200}^{\mathrm{int}}`$ is tightly correlated with mass and mass-weighted temperature in all runs. As we shall see, however, these correlations show significantly less scatter than those relating $`Y`$ with X-ray emission-weighted properties, because the latter is sensitive to substructure in the dense core gas. The best power-law fits to the $`Y_{200}^{\mathrm{int}}M_{200}`$ relation in each simulation are

* Non-radiative run:

$$Y_{200}^{\mathrm{int}}=1.22\times 10^6\left(\frac{M_{200}}{10^{14}h^1\mathrm{M}_{}}\right)^{1.69}\left(h^1\mathrm{Mpc}\right)^2,$$
(11)
* Cooling run:

$$Y_{200}^{\mathrm{int}}=0.97\times 10^6\left(\frac{M_{200}}{10^{14}h^1\mathrm{M}_{}}\right)^{1.79}\left(h^1\mathrm{Mpc}\right)^2,$$
(12)
* Preheating run:

$$Y_{200}^{\mathrm{int}}=1.12\times 10^6\left(\frac{M_{200}}{10^{14}h^1\mathrm{M}_{}}\right)^{1.93}\left(h^1\mathrm{Mpc}\right)^2.$$
(13)

As expected, the slope of the fit to the Non-radiative clusters is close to 5/3, the self-similar value given in Eq. (6). The slope is steeper for the Cooling and Preheating simulations however, primarily due to lower $`Y`$ values for low-mass systems. This difference is due to the effects of cooling and heating on both the gas fraction and temperature of the clusters. Cooling removes low-entropy (core) gas from the halo, causing higher-entropy material to flow in to replace it. This results in a lower gas fraction although the remaining gas is hotter (Pearce et al. 2000; Bryan 2000; Muanwong et al. 2001). For low-mass clusters ($`M_{200}<10^{14}h^1`$M), the net effect of cooling on $`Y`$ is dominated by the decrease in hot gas fraction, shown in Fig. 2, around a factor of two lower than the Non-radiative average at $`M_{200}=10^{13}h^1`$M. For the highest mass clusters this difference decreases to around $`30`$ per cent; combined with the increase in temperature of the gas, the clusters have very similar $`Y`$ values in the two simulations.

In the Preheating case, the energy injection at high redshift heats the gas and expels some of it from the cluster. Hence its effect on the cluster properties is similar to the cooling model but the fate of the gas is different (Muanwong et al. 2002). Fig. 2 shows that the gas fraction in Preheating clusters increases more rapidly with mass than for objects in the Cooling simulation (please refer to Section 3.1 and Fig. 3 of Muanwong et al. 2002 for a detailed discussion of the hot, cold and total baryonic mass fractions in clusters for the same set of simulations). In low-mass systems the heating energy was sufficient to significantly decrease $`f_{\mathrm{gas}}`$ (and hence $`Y`$), but left gas fractions almost unchanged in the highest-mass clusters, where the heating resulted only in a modest increase in temperature. Note, therefore, that the best-fit relation is only a good description of the low-mass clusters in the Preheating simulation; systems more massive than $`10^{14}h^1M_{}`$ are adequately described by the Non-radiative relation. This also explains why White et al. (2002) find good agreement with the self-similar $`Y_{\mathrm{int}}M`$ scaling, as they only investigated high-mass clusters ($`M_{200}5\times 10^{14}h^1\mathrm{M}_{}`$), where a high fraction of the wind material would have been retained.

The effects of cooling and preheating on the cluster properties are more readily apparent in the relation between $`Y`$ and mass-weighted temperature. For a cluster of fixed gas and total mass, $`YT_{\mathrm{mw}}`$, and so heating the gas causes a shift to the right of the steeper self-similar relation ($`YT^{5/2}`$). Additionally, the decrease in gas fraction lowers $`Y`$, separating the two relations further. The best-fit $`Y_{200}^{\mathrm{int}}T_{\mathrm{mw}}`$ relations are

* Non-radiative run:

$$Y_{200}^{\mathrm{int}}=5.59\times 10^7\left(\frac{k_\mathrm{B}T_{\mathrm{mw}}}{1\mathrm{keV}}\right)^{2.63}\left(h^1\mathrm{Mpc}\right)^2$$
(14)
* Cooling run:

$$Y_{200}^{\mathrm{int}}=2.51\times 10^7\left(\frac{k_\mathrm{B}T_{\mathrm{mw}}}{1\mathrm{keV}}\right)^{2.90}\left(h^1\mathrm{Mpc}\right)^2$$
(15)
* Preheating run:

$$Y_{200}^{\mathrm{int}}=2.23\times 10^7\left(\frac{k_\mathrm{B}T_{\mathrm{mw}}}{1\mathrm{keV}}\right)^{3.03}\left(h^1\mathrm{Mpc}\right)^2.$$
(16)

### 4.2 The $`YT_\mathrm{X}`$ relation

Of more practical interest is to investigate how $`Y`$ correlates with X-ray observables, particularly luminosity and emission-weighted temperature.

Fig. 3 illustrates the relationship between $`Y_{200}^{\mathrm{int}}`$ and $`T_\mathrm{X}`$ for the Cooling (top panel) and Preheating (bottom panel) simulations. Circles represent clusters for which all hot gas particles within $`R_{200}`$ are considered, whereas crosses are for the same clusters but excluding particles within the cooling radius in the computation of the temperature. We define the cooling radius as the radius within which the gas has a mean cooling time of 6 Gyr. (This procedure is an attempt to provide a ‘cooling flow’ correction to the estimated values of $`T_\mathrm{X}`$.)

Comparing these results with Fig. 1, we see that the scatter is larger for the correlation with $`T_\mathrm{X}`$ than with $`T_{\mathrm{mw}}`$. This is because $`T_\mathrm{X}`$ is weighted by density and temperature and is therefore more sensitive to substructure, particularly from the centres of clusters. This is also evident from the fact that the scatter is smaller when the gas within the cooling radius is excluded (note that this also increases $`T_\mathrm{X}`$ in high-mass systems, which slightly flattens the $`YT_\mathrm{X}`$ relation).

Table 1 lists best-fit coefficients, $`\alpha _\mathrm{T}`$ and $`A_\mathrm{T}`$, to the power-law, $`Y_{200}^{\mathrm{int}}=(A_\mathrm{T}/10^7)(k_\mathrm{B}T_\mathrm{X}/1\mathrm{k}\mathrm{e}\mathrm{V})^{\alpha _\mathrm{T}}`$ for the Cooling and Preheating runs. (We also present results from the Non-Radiative run, using a bolometric cooling function, and this agrees well with the self-similar scaling $`YT^{5/2}`$.) For a given run, the normalization of the $`Y_{200}^{\mathrm{int}}T_\mathrm{X}`$ relation is lower and slope steeper than for the $`Y_{200}^{\mathrm{int}}T_{\mathrm{mw}}`$ relation. The heavier weighting of the denser gas exacerbates the effects of cooling/heating on the $`YT`$ relation discussed in the previous section.

### 4.3 The $`YL_\mathrm{X}`$ relation

We end this section by reporting on the $`z=0`$ correlation between $`Y^{\mathrm{int}}`$ and bolometric X-ray luminosity, $`L_\mathrm{X}`$, as inferred from the soft-band emission. Fig. 4 shows the $`Y_{200}^{\mathrm{int}}L_\mathrm{X}`$ relation from the Cooling (top panel) and Preheating (bottom panel) simulations. Again, we use circles to represent values including all hot gas particles, whereas crosses represent results when excluding gas from within the cooling radius. The best-fit parameters $`\alpha _\mathrm{L}`$ and $`A_\mathrm{L}`$, defined by $`Y_{200}^{\mathrm{int}}=(A_\mathrm{L}/10^7)(L_\mathrm{X}/L_{43})^{\alpha _\mathrm{L}}`$, where $`L_{43}=10^{43}h^2\mathrm{erg}\mathrm{s}^1`$, are listed in Table 1.

The Cooling and Preheating relations are significantly shallower than the predicted self-similar scaling, $`YL^{5/4}`$, which is well reproduced by the Non-radiative simulation (where $`\alpha _\mathrm{L}=1.26`$). Again, this is due to the increase in entropy of the gas, particularly in low-mass clusters, which decreases their luminosity. Omitting gas from within the cooling radius slightly increases the slope of the $`Y^{\mathrm{int}}L_\mathrm{X}`$ relation and also reduces the scatter.

## 5 Evolution of scaling relations

We now investigate the evolution of the scaling relations with redshift. In each of the following subsections, we factor out the dependence expected from self-similar evolution. Note that the self-similar evolution has no dependence on $`\mathrm{\Omega }_0`$ because we define the overdensity with respect to the comoving critical density.

### 5.1 Evolution of the $`YM`$ relation

Fig. 5 shows the $`Y_{200}M_{200}`$ relation for clusters in the Cooling simulations, for a range of redshifts between $`z0`$ and $`z1`$. Here we use $`Y_{200}=Y_{200}^{\mathrm{int}}d_\mathrm{A}^2`$, demonstrating the dominant effect of $`d_\mathrm{A}`$ on $`Y`$ at these redshifts. Note that the number of clusters satisfying our mass selection criteria, $`M_{200}1.18\times 10^{13}h^1\mathrm{M}_{}`$, decreases with increasing redshift.

A noticeable feature of Fig. 5 is the apparent independence of the slope of the $`Y_{200}M_{200}`$ relation with redshift. To quantify the effects of cooling and heating on the evolution of the slope, we fit a power-law

$$\frac{Y_{200}^{\mathrm{int}}}{1+z}=10^{\beta _\mathrm{M}(z)}(M_{200}/10^{14}h^1M_{})^{\alpha _\mathrm{M}(z)},$$
(17)
to the distribution of clusters at each redshift, and plot the resulting best-fit values of $`\alpha _\mathrm{M}`$ against $`z`$ in the upper panel of Fig. 6. (We do not plot results for redshifts beyond 2 as the number of clusters below our mass threshold becomes too small for a reliable fit to be produced.) Crosses, circles and diamonds represent best-fit values, and dotted lines denote 1-$`\sigma `$ uncertainties (estimated using the method given by Press et al. 1992) from the Non-radiative, Cooling and Preheating simulations respectively. Indeed, the results are consistent with an unevolving slope for all three models out to $`z=2`$.

The lower panel of Fig. 6 shows the evolution of the normalization, obtained by fitting the same power-law function as before, but with a fixed slope taken to be equal to the mean value of $`\alpha _\mathrm{M}(z)`$ within the displayed range of $`z`$ (values are given in Table 2). From this figure, we see that the evolution of the $`Y_{200}^{\mathrm{int}}M_{200}`$ normalization is remarkably self-similar for all runs. Fitting a power-law, $`10^{\beta _\mathrm{M}(z)}(1+z)^{\gamma _\mathrm{M}}`$, to the derived normalizations (within the displayed range of $`z`$), we find that the largest deviation from self-similar evolution of the normalization is only $`\gamma _\mathrm{M}=0.14`$, for the Cooling simulation (Table 2).

It is therefore clear that the evolution of the $`Y_{\mathrm{int}}M_{200}`$ relation is not a sensitive probe of underlying non-gravitational physics, as it is a measure of the large-scale evolution of clusters. As discussed in McCarthy et al. (2003a), the effects of cooling and preheating are much more apparent if the SZ flux density is evaluated within a smaller radius, where the average gas density is higher. As we shall see, such effects are apparent when studying evolution of the X-ray properties of the clusters.

### 5.2 Evolution of the $`YT_\mathrm{X}`$ relation

Fig. 7 shows the redshift dependence of the slope and normalization of the $`YT_\mathrm{X}`$ relation for the Cooling and Preheating simulations, using cooling-flow corrected temperatures. As in the previous section, we first fit the cluster distribution at each redshift with a power-law

$$\frac{Y_{200}^{\mathrm{int}}}{(1+z)^{3/2}}=10^{\beta _\mathrm{T}(z)}(T_\mathrm{X}/\mathrm{keV})^{\alpha _\mathrm{T}(z)},$$
(18)
to determine the evolution of the slope. While this shows significant statistical oscillations (reflecting the higher sensitivity of the emission-weighted temperature to the thermal history of the gas), there is no evidence of evolution with redshift beyond that predicted by self-similarity. The slope is then fixed at the mean value, $`\alpha _\mathrm{T}`$, before determining the redshift dependence of the normalization, $`\beta _\mathrm{T}(z)`$. Unlike the $`Y^{\mathrm{int}}`$$`M`$ relation, there is a definite evolutionary trend that is particularly significant in the Cooling case. This is an indication that evolution is predominantly in the X-ray properties rather then in the SZ signal, which is much less sensitive to the details of the gas distribution in the central regions of clusters. Fitting the normalizations with a power-law of the form $`10^{\beta _\mathrm{T}(z)}(1+z)^{\gamma _\mathrm{T}}`$, we find $`\gamma _\mathrm{T}=0.73`$ and $`\gamma _\mathrm{T}=0.33`$ for the Cooling and Preheating simulations, respectively. In both cases the normalization of the $`T_\mathrm{X}`$$`M`$ relation increases with time compared to the expected self-similar evolution.

### 5.3 Evolution of the $`YL_\mathrm{X}`$ relation

We end by studying the evolution of the $`YL_\mathrm{X}`$ relation, shown in Fig. 8. Our power-law fitting function was

$$\frac{Y_{200}^{\mathrm{int}}}{(1+z)^{27/8}}=10^{\beta _\mathrm{L}(z)}(L_\mathrm{X}/L_{43})^{\alpha _\mathrm{L}(z)},$$
(19)
using cooling-flow corrected values for $`L_\mathrm{X}`$. Again, the slopes are constant with redshift to good approximation, and there is additional evolution in the normalization. In this case a single power-law is not so good a fit over the complete redshift range, and we restrict the fit to $`0<z<1`$ which is where observations are likely to be made. Fitting a power-law $`10^{\beta _\mathrm{L}(z)}(1+z)^{\gamma _\mathrm{L}}`$ yields $`\gamma _\mathrm{L}1.3`$ & $`1.6`$ in the Cooling and Preheating cases respectively (Table 2) in this redshift range, with the relation flattening at redshifts beyond one.

## 6 Conclusions

In this paper, we investigated the relationship between SZ and X-ray properties of clusters drawn from three simulations: a Non-radiative simulation, where the gas could only heat adiabatically and through shocks; a Cooling simulation, where the gas could also cool radiatively and a Preheating simulation, where we additionally heated the gas by 1.5 keV per particle at $`z=4`$. We focussed on the relations between the thermal SZ signal $`Y_{\mathrm{int}}`$ and cluster mass, temperature and X-ray luminosity, both at $`z=0`$ and the evolution with redshift. Our main conclusions are as follows:

* $`Y_{\mathrm{int}}`$ is tightly correlated with both mass and mass-weighted temperature. Scaling relations of $`Y_{\mathrm{int}}`$ with X-ray quantities show significantly more scatter however, due to the sensitivity of the latter (which are emission-weighted) to the gas distribution in the cores of clusters. Excluding gas from within the cooling radius of each cluster reduces the scatter.
* The Non-radiative simulation scaling relations closely match those predicted by the self-similar model, with $`Y_{\mathrm{int}}T_\mathrm{X}^{5/2}`$ and $`Y_{\mathrm{int}}L_\mathrm{X}^{5/4}`$.
* Cooling decreases $`Y_{\mathrm{int}}`$ for a given cluster mass, due to the reduction in the hot gas fraction (which decreases $`Y_{\mathrm{int}}`$) being more effective than the increase in temperature of the remaining gas (which increases $`Y_{\mathrm{int}}`$). Preheating the gas has a similar effect except in the most massive clusters, where the heating energy was not sufficient to expel much gas from the halo. At redshift zero, the resulting $`Y_{\mathrm{int}}`$$`T_\mathrm{X}`$ scaling relation is steeper than the self-similar relation, $`Y_{\mathrm{int}}T_\mathrm{X}^{3.03.4}`$, and the $`Y_{\mathrm{int}}`$$`L_\mathrm{X}`$ relation flatter, $`Y_{\mathrm{int}}L_\mathrm{X}^{0.91.1}`$.
* There is no evidence of evolution in the slopes of any of the scaling relations. The normalization of the $`Y_{200}M_{200}`$ relation is also consistent with self-similar evolution. However, the relations between $`Y_{\mathrm{int}}`$ and the cluster X-ray properties do show significant evolution in normalization relative to that expected from self-similarity, indicating that evolution is predominantly in the X-ray properties rather than the SZ properties. The effects of excess entropy reduces the negative evolution of the $`Y_{\mathrm{int}}`$$`T_\mathrm{X}`$ and $`Y_{\mathrm{int}}`$$`L_\mathrm{X}`$ with redshift, particularly at $`z<1`$.

## Acknowledgments

A.C.S. was supported by the EU CMBNET network, S.T.K. by PPARC, and A.R.L. in part by the Leverhulme Trust. The simulations used in this paper were carried out on the Cray-T3E (RIP) at the EPCC as part of the Virgo Consortium programme of investigations into the formation of structure in the Universe. A.C.S. acknowledges the use of the computer facilities at CALMIP-CICT, France, and visits to Sussex supported by PPARC.


format: text

text:
# Pulse Propagation in Chains with Nonlinear Interactions

## Abstract

Pulse propagation in nonlinear arrays continues to be of interest because it provides a possible mechanism for energy transfer with little dispersion. Here we show that common measures of pulse dispersion might be misleading; in strongly anharmonic systems they tend to reflect a succession of extremely narrow pulses traveling at decreasing velocities rather than the actual width of a single pulse. We present analytic estimates for the fraction of the initial energy that travels in the leading pulses. We also provide analytic predictions for the leading pulse velocity in a Fermi-Pasta-Ulam $`\beta `$-chain.

The stability of localized energy, e.g. in the form of breathers, in translationally invariant nonlinear arrays, and the way in which localized energy packets can be transported in these arrays, has been a topic of interest for several decades, and continues to be of great interest for a number of reasons. One is that many of the ideas on the subject have recently and increasingly been confirmed experimentally. Another is the possible importance of the subject in the transport of energy in biological systems. A third is the ever increasing numerical capability that allows simulations of larger systems over longer times. A recent focus issue of the journal *Chaos* chaos  contains some of the most current contributions and reviews of the subject, and covers the three topics just mentioned.

While the advances of the past few years are exciting and enormously instructive, the analytic understanding of these phenomena has been made difficult by the fact that the systems are nonlinear. Many of the available results (including those obtained in our group) are numerical, and it is sometimes difficult and even misleading to generalize from these results (for reviews of the subject preceding the special issue noted above, see flach ; aubry  and references therein). Our contribution in this paper is an analytic understanding of results previously obtained only numerically.

A typical set of questions that one can pose is the following: Suppose that a single unit in a nonlinear array is given an initial velocity. How will this velocity/energy propagate through the array? Will some or all of the energy remain localized, or will it spread? If a localized moving pulse does develop, at what velocity will it propagate? These are some of the signal propagation issues that we address analytically.

We focus on the one-dimensional Fermi-Pasta-Ulam (FPU)-type problem for unit mass particles described by the Hamiltonian

$$H=\underset{i}{}\frac{\dot{x}_i^2}{2}+\underset{i}{}V(x_ix_{i1}),$$
(1)
where $`x_i`$ is the displacement of particle $`i`$ from its equilibrium position, and $`V(z)`$ is the potential

$$V_n(z)=\frac{k^{}}{n}\underset{i}{}|z|^n+\frac{k}{2}\underset{i}{}z^2.$$
(2)
For the FPU $`\beta `$-problem $`n=4`$, but we retain $`n`$ as a general power because a number of theories and simulations deal with other values of $`n`$, and portions of our analysis do as well. The parameters $`k`$ and $`k^{}`$ are the harmonic and anharmonic force constants, respectively. The variables and the time can be scaled so that the only distinct cases of this problem are $`k=0`$ (purely anharmonic chain), $`k^{}=0`$ (purely harmonic chain), and $`k,k^{}0`$ (“mixed” chain; $`k=k^{}=1`$ is a convenient choice). The control parameter is then the initial velocity. The equation of motion for the $`i^{th}`$ particle in a mixed chain then is

$$\begin{array}{cc}\hfill \ddot{x}_i& =\left|x_{i+1}x_i\right|^{n1}\mathrm{sgn}(x_{i+1}x_i)+(x_{i+1}x_i)\hfill \\ & +\left|x_ix_{i1}\right|^{n1}\mathrm{sgn}(x_ix_{i1})+(x_ix_{i1}),\hfill \end{array}$$
(3)
where $`\mathrm{sgn}(x)=\pm 1`$ for $`x0`$. Initially all particles are at rest in their equilibrium positions except for one particle (far from any boundaries) that has initial velocity $`v_0`$. We take the chain to be sufficiently long and the boundaries sufficiently far from the initial excitation that their precise nature does not matter for our analysis.

In sarmiento , Sarmiento et al. analyze the pulse evolution in terms of the mean distance from the initial site (“pulse position”) and its dispersion (“pulse width”),

$$x=\frac{_i\left|i\right|E_i}{_iE_i},\sigma ^2=\frac{_ii^2E_i}{_iE_i}x^2.$$
(4)
The local energy is defined as

$$E_i=\frac{\dot{x}_i^2}{2}+\frac{1}{2}V(x_{i+1},x_i)+\frac{1}{2}V(x_i,x_{i1}).$$
(5)
We begin by considering the spreading of the initial pulse. The width of the pulse as time proceeds is often invoked as a measure of the ability of the nonlinearity to keep the energy localized. It is well known that in a harmonic lattice an initial pulse spreads even as it moves. In a mixed chain, one expects less spreading for a higher initial velocity $`v_0`$ since a more energetic pulse samples the more anharmonic portions of the potential. Indeed, highly localized breathers have been shown to be exact solutions for the purely anharmonic chain in the limit $`n\mathrm{}`$ page . For a quartic anharmonicity ($`n=4`$) the contribution of the harmonic and anharmonic contributions to the potential energy are equal at the maximum displacement associated with kinetic energy $`v_0^2/2=4`$. Therefore when $`v_0\sqrt{8}`$ ($`v_0\sqrt{8}`$) the dominant contribution to the potential energy of the pulse is the harmonic (anharmonic) portion.

In Fig. 1 we show the pulse width as a function of pulse position for $`n=4`$ in two cases. In one, $`v_0=1`$ so that the harmonic portion of the potential is strongly sampled by the excitation. In the other, $`v_0\mathrm{}`$, i.e., the potential is essentially a purely quartic potential. The pulse in the purely anharmonic potential is more localized after traveling a given mean distance than is the pulse in the mixed system sarmiento . The inset shows the mean pulse position as a function of time (the pulses move at a constant speed).

We note that the lower energy pulse moves more rapidly than the higher energy pulse. This result may appear contradictory with those obtained earlier in sarmiento  where, *for a given initial pulse energy*, pulse velocities in purely harmonic and purely anharmonic systems were compared. A pulse in a purely harmonic system moves at a speed that is independent of the initial velocity $`v_0`$, while the pulse speed in a purely quartic system increases with increasing $`v_0`$ sarmiento . Therefore, for sufficiently high initial pulse velocity, a pulse in a purely harmonic system moves more slowly than in a purely anharmonic chain sarmiento . In a mixed chain, on the other hand, the pulse speed as a function of $`v_0`$ is bounded *below* by the higher of the two (purely harmonic and purely quartic), approaching the purely harmonic behavior at low $`v_0`$ and the purely quartic behavior at large $`v_0`$.

Alternatively but equivalently, one can say that the velocity of a pulse in a purely harmonic chain is independent of its amplitude while in an anharmonic chain the pulse speed decreases with decreasing amplitude. This observation has important implications in the understanding of the way in which the pulse width increases as the pulse travels away from the site of origin. To see this, consider a low-velocity pulse launched in a mixed chain. The harmonic portion of the potential here dominates the evolution of the pulse, as shown in Fig. 2(a). Two symmetric fronts travel away from the origin carrying part of the energy and spreading. The remaining energy is progressively distributed among the particles between the pulses. Altogether one observes a gradual broadening of the two pulses, exactly as one imagines a pulse broadening to occur.

However, for a mixed potential with high initial pulse velocity the situation is quite different, as can be seen in Fig. 2(b). Now a portion of the energy travels symmetrically away from the center in extremely localized pulses that in fact *remain highly localized*. The remaining energy is successively “launched” from the origin in the form of secondary pulses of smaller and smaller amplitude which travel more and more slowly (the low initial velocity or nearly harmonic case can be thought of in these terms as well, but the secondary pulses travel at essentially the same speed as the primary pulses). Hence, there is a series of narrow pulses of decreasing amplitude that are getting further apart from one another, giving rise to the apparent “dispersion”. Moreover, the maxima of the pulses oscillate; when the energy is more concentrated in one particle the maximum is higher than when it is shared between two or more.

The extreme localization of the pulses suggests that a two- or three-particle approximation may capture the essence of the physics of the problem, and it is this feature that we use to arrive at analytic results. The point to stress here is that *the second moment* $`\sigma ^2`$ *is in fact a deceptive measure of the dispersion of the pulse* except in an essentially harmonic system.

We begin by estimating the energy in the primary pulses and, from this, the pulse velocity as determined by the pulse energy. We assume that the only effects of the restoring forces are to split part of the energy into two pulses, and to keep the remainder of the energy at the origin, from where it will create the secondary pulses. Therefore we need to calculate how much energy is transmitted from the particle at $`i=0`$ to the particles at $`i=\pm 1`$, and how quickly it is transmitted.

The first step is achieved by considering a three particle system and neglecting the rest of the chain. This approximation presupposes that the potential is sufficiently steep that the particles at $`i=\pm 2`$ barely move before the particles at $`i=\pm 1`$ have acquired their full velocity. Obviously, this is not strictly true for finite $`n`$; however, we will show that it is a very good approximation, especially for high pulse velocities.

Initially, the three particles are in their equilibrium positions and all of the energy $`E_0=v_0^2/2`$ is concentrated in the middle particle, $`i=0`$. Some of the energy is transfered to the neighbors as the springs compress and stretch, and at some time later the energy of the three particles is once again all kinetic. The symmetry of the system requires that this occur when the particles at $`i=\pm 1`$ acquire their maximum velocity, $`u`$ (equal by symmetry), and the one at $`i=0`$ its minimum velocity, $`u_0`$. The three-particle system oscillates back and forth between these two configurations, but we are only interested in this first portion of the cycle. Energy and momentum conservation lead to

$$u=\frac{2}{3}v_0,u_0=\frac{1}{3}v_0.$$
(6)
Each primary pulse therefore carries away an energy

$$E_p=\frac{1}{2}u^2=\frac{2}{9}v_0^2=\frac{4}{9}E_0.$$
(7)
Figure 3 shows simulation results for the primary pulse energy as a function of the power $`n`$ of the potential for two cases. The circles correspond to an initial pulse velocity $`v_0\mathrm{}`$ (or, alternatively, a purely quartic potential with any $`v_0`$). The squares correspond to an initial pulse velocity $`v_0=10`$. Both lie in the regime where the dynamics is dominated by the anharmonicity, $`v_0\sqrt{8}`$, but our theory is expected to improve with increasing $`v_0`$. The asymptotic value for the high $`v_0`$ case is $`E_p/E_0=4/9`$, exactly as predicted. The agreement is very good even for $`n=4`$, where our prediction is already within a few percent of the correct value. For the lower initial pulse velocity the energy ratio is asymptotically only about 7% smaller than predicted, a reflection of the fact that as the harmonic component becomes more important the pulse occupies more than three sites.

Next we turn our attention to the calculation of the velocity of the primary pulse when $`n=4`$ (the FPU case). For this calculation we simplify our model even further and consider only a two-particles system, one of which has the initial velocity $`u=2v_0/3`$. We then calculate the time, $`T(v_0)`$, for the second particle to acquire the same velocity as the first, i.e., the time at which the velocities of the two particles are equal. We maintain that this is the time necessary for the primary pulse to travel from one particle to the next. To calculate this time it is not necessary to actually integrate the equations of motion. Defining $`z=x_1x_2`$, we have for $`n=4`$, $`\ddot{z}=2z^32z`$. This equation of motion is derived from a potential $`V(z)=z^4/2+z^2`$, and the initial conditions are $`z(0)=0`$ and $`\dot{z}(0)=u`$. From energy conservation we have that the final energy (kinetic plus potential) is equal to the initial energy (all kinetic),

$$\frac{1}{2}\dot{z}^2+V(z)=\frac{1}{2}u^2,$$
(8)
from which it follows that $`\dot{z}=\sqrt{u^22z^2z^4}`$. Furthermore, the two particles have the same velocity when $`\dot{z}=0`$. It immediately follows that the relative displacement at that moment is given by $`z_m=\sqrt{\sqrt{1+u^2}1}`$. Now we can integrate $`\dot{z}`$,

$$T(v_0)=\frac{1}{z_m}_0^1\frac{\mathrm{d}\eta }{\sqrt{\left(1\eta ^2\right)\left(\mathrm{coth}^2(\varphi /2)+\eta ^2\right)}},$$
(9)
where we have introduced the variables $`\eta z/z_m`$ and $`\varphi =\mathrm{sinh}^1(u)`$. This integral can be done analytically grad :

$$T(v_0)=\frac{1}{\sqrt{2}\mathrm{cosh}^{1/2}(\varphi )}K\left(\frac{\left(\mathrm{cosh}(\varphi )1\right)^{1/2}}{\sqrt{2}\mathrm{cosh}^{1/2}(\varphi )}\right),$$
(10)
where $`K(x)`$ is the complete elliptic integral of the first kind grad . The pulse velocity is just the inverse of this time:

$$C(v_0)=T^1(v_0).$$
(11)
For $`v_01`$ (strongly anharmonic potential) one finds from dimensional analysis ourchain  that $`C(v_0)v_0^{1/2}`$.

In Fig. 4, we compare the results of our approximation with the numerical simulation of the full chain as a function of $`v_0`$. The circles are the simulation results, and the broken line is our two-particle approximation, Eq. (11). The agreement is clearly excellent for initial pulse velocities above $`v_0>\sqrt{8}`$, the value that we offered as a limit for the validity of this approach.

In this communication we have shown that traditional measures of pulse propagation in arrays with nonlinear interactions may be misleading. We discussed systems with interactions that have both harmonic and anharmonic contributions, and argued that the behavior of a pulse launched by imparting an initial velocity $`v_0`$ to a particle at one site in such a system depends strongly on $`v_0`$. For low velocities the harmonic portions of the potential are primarily sampled, and the pulse behaves as it would in a harmonic system. For high velocities the anharmonic portions of the potential dominate the behavior, and the pulse propagates as it would in a purely anharmonic chain. It is in the anharmonic regime that one must view traditional measures of pulse propagation with some caution. In particular, we have shown that in the anharmonic regime the usual second moment “pulse width” is not a measure of the way a single pulse spreads, but rather, of the span covered by a series of very narrow pulses of decreasing velocity. In a statistical measure this appears as a growing second moment. We have also presented analytic estimates for the energy and velocity of the leading pulse, and have shown by comparison with numerical simulations that our estimates are extremely accurate in the anharmonic regime.

A number of extensions of our approach are possible, albeit with some analytic complications. For example, the approach can be extended to the FPU $`\alpha `$-model that includes cubic as well as quartic interactions. Some of the polynomial solutions and integrations that we have carried out analytically might then have to be done numerically. The same is true were one to include dissipation. The model can be extended to include a local harmonic potential, and also to higher dimensions. We are currently exploring these generalizations.

This work was supported by the Engineering Research Program of the Office of Basic Energy Sciences at the U. S. Department of Energy under Grant No. DE-FG03-86ER13606.


text:
# Exact solitary wave solutions of the nonlinear Schrödinger equation with a source

## Abstract

We use a fractional transformation to connect the traveling wave solutions of the nonlinear Schrödinger equation (NLSE), phase-locked with a source, to the elliptic functions satisfying, $`f^{\prime \prime }\pm af\pm \lambda f^3=0`$. The solutions are necessarily of the rational form, containing both trigonometric and hyperbolic types as special cases. Bright, and dark solitons, as also singular solitons, are obtained in suitable range of parameter values.

Much attention has been paid to the study of the externally driven NLSE, after the seminal work of Kaup and Newell kaup . This equation features prominently in the problem of optical pulse propagation in asymmetric, twin-core optical fibers gil ; synder ; boris ; currently an area of active research. Of the several applications of an externally driven NLSE, perhaps the most important ones are to long Josephson junctions sam , charge density waves newell , and the plasmas driven by rf fields bekki . The phenomenon of autoresonance land , indicating a continuous phase locking between the solutions of NLSE and the driving field has been found to be a key characteristic of this system. In the presence of damping, this dynamical system exhibits rich structure including bifurcation. This is evident from analyses, around a constant background, as well as numerical investigations smirnov ; barashenkov ; nistazakis . Although the NLSE is a well-studied integrable system agrawal , no exact solutions have so far been found for the NLSE with a source, to the best of the authors’ knowledge. All the above inferences have been drawn through perturbations around solitons and numerical techniques.

In this Letter, we map exactly, the traveling wave solutions of the NLSE phase-locked with a source, to the elliptic functions, through a fractional transformation (FT). It was found that the solutions are necessarily of the rational type, with both the numerator and denominator containing terms quadratic in elliptic functions, in addition to having constant terms. It is well-known that the solitary wave solutions of the NLSE whitham ; novikov  are cnoidal waves, which contain the localized soliton solutions in the limit, when the modulus parameter equals one abramowitz . Hence, the solutions found here, for the NLSE with a source, are nonperturbative in nature. We find both bright and dark solitons as also singular ones. Solitons and solitary pulses show distinct behavior. In the case, when the source and the solutions are not phase matched, perturbation around these solutions may provide a better starting point.

For nonlinear equations, a number of transformations are well-known in the literature, which map the solutions of a given equation to the other drazin ; das . The familiar example being the Miura transformation miura , which maps the solutions of the modified KdV to those of the KdV equation. To find static and propagating solutions, appropriate transformations have also been cleverly employed, to connect the nonlinear equations to the ones satisfied by the elliptic functions: $`f^{\prime \prime }\pm af\pm \lambda f^3=0`$. Here and henceforth, prime denotes derivative with respect to the argument of the function. Solitons and solitary wave solutions of KdV, NLSE, and sine-Gordon etc., can be easily obtained in terms of the elliptic functions in this manner.

We consider the NLSE with an external traveling wave source:

$$i\frac{\psi }{t}+\frac{^2\psi }{x^2}+g\psi ^2\psi +\mu \psi =ke^{i[\chi (\xi )\omega t]},$$
(1)
where $`g`$, $`\mu `$, and $`k`$ are real and $`\xi =\alpha (xvt)`$. The traveling wave solution, phase-locked with the external source is taken to be $`\psi (x,t)=e^{i[\chi (\xi )\omega t]}\rho (\xi ).`$ Separating the real and imaginary parts of Eq. (1), and integrating the imaginary part one gets,

$$\chi ^{}=\frac{v}{2\alpha }+\frac{c}{\alpha \rho ^2},$$
(2)
where $`c`$ is the integration constant. In order that, the external phase is independent of $`\psi `$, we put $`c=0`$ to obtain,

$$\alpha ^2\rho ^{\prime \prime }+g\rho ^3ϵ\rho k=0,$$
(3)
where $`ϵ=\omega \frac{v^2}{4}\mu `$. For illustrating our procedure, we assume that $`ϵ>0`$, and $`g>0`$, the other cases can be studied analogously. After a suitable scaling of the field variable, the above equation can be cast into the convenient form,

$$\rho ^{\prime \prime }+\rho ^3\rho k=0.$$
(4)
We find, after a straightforward but lengthy algebra that, the following FT,

$$\rho (\xi )=\frac{A+Bf^\delta (\xi ,m)}{1+Df^\delta (\xi ,m)},$$
(5)
for $`ABD0`$, maps the solutions of Eq. (4) to the elliptic functions $`f(\xi )`$, provided $`\delta =2`$. $`f(\xi )`$ satisfies $`f^{\prime \prime }=ff^3`$, with a conserved quantity $`E_0=f^2/2+(1/4)f^4f^2/2`$. The above equation has four bounded periodic solutions which are Jacobi elliptic functions: $`\mathrm{cn}(\xi ,m)`$, $`\mathrm{sd}(\xi ,m)`$, $`\mathrm{dn}(\xi ,m)`$, and $`\mathrm{nd}(\xi ,m)`$, where $`m`$ is the modulus parameter. Apart from finite energy oscillatory solutions, the above has localized soliton solutions. It should be noted that for the attractive case $`(g<0)`$, the bounded solutions are $`\mathrm{sn}(\xi ,m)`$, and $`\mathrm{cd}(\xi ,m)`$. For $`AD=B`$, one only finds a constant amplitude solution.

Before elaborating on specific cases, a number of interesting features emerging from the above mapping is worth mentioning. First of all, the nontrivial solutions are necessarily of the rational type, i.e., $`D0`$, in the presence of the source. Secondly, $`E_0=0`$, and $`E_00`$ cases, show characteristically different behavior. For example, for $`E_00`$ case, one can have solutions with $`A=0`$, and $`B0`$, of the type

$$\rho (\xi )=\frac{(k/4E_0)f^2}{1+(1/8E_0)f^2}.$$
(6)
However, $`A0`$, and $`B=0`$ is not allowed. For $`E_0=0`$ case, one can have for $`B=0`$, and $`A0`$, a singular solution of the following form,

$$\rho (\xi )=\frac{2k}{1f^2}.$$
(7)
As we will see explicitly below, nonsingular solutions are also possible. Like the former case, here $`B0`$, and $`A=0`$ is not allowed in the presence of the source.

For explicitness, we consider the unscaled equation containing all the parameters and illustrate below various type of solutions, taking $`f=cn(\xi ,m)`$. Other cases can be similarly worked out. The consistency conditions are given by,

$`Aϵ2\alpha ^2(ADB)(1m)+gA^3k`$ $`=0,`$ (8)
$`2ϵAD+ϵB+6\alpha ^2(ADB)D(1m)`$
$`4\alpha ^2(ADB)(2m1)+3gA^2B3kD`$ $`=`$ $`0,`$ (9)
$`AϵD^2+2ϵBD+4\alpha ^2(ADB)D(2m1)+`$
$`6\alpha ^2(ADB)m+3gAB^23kD^2`$ $`=`$ $`0,`$ (10)
$`ϵBD^22\alpha ^2(ADB)Dm+gB^3kD^3`$ $`=`$ $`0.`$ (11)
The above equations clearly indicate that the solutions, for $`m=1`$, $`m=0`$ and other values of $`m`$, have distinct properties. For example, when $`m=1`$, $`A`$ is obtained as the solution of the cubic equation \[Eq. (8)\], containing the source strength $`k`$. Similarly, for $`m=0`$, either $`B`$, or $`D`$ appears as the solution of Eq. (11). As noted earlier, when $`D=0`$, $`B`$ also equals zero, indicating only a constant solution $`(A0)`$. One needs to be careful in choosing the real solutions of the above cubic equations, for suitable range of the parameter values. Although a wide class of solutions are allowed; for brevity, we only outline a few of the interesting solutions and their properties in a restricted range of equation parameters.

We start with the localized solitons, by taking $`m=1`$ $`(\mathrm{cn}(\xi ,1)=\mathrm{sech}(\xi ))`$, with the parameter values $`A=1`$, $`B=\delta `$, and $`D=1\delta `$ we obtain

$$\rho (\xi )=\frac{1+\delta \mathrm{sech}^2(\xi )}{1+\mathrm{\Gamma }\mathrm{sech}^2(\xi )};$$
(12)
here the amplitude, width, and velocity are related as, $`\delta =(Q+\sqrt{Q^24PR})/2Q`$, with $`Q=ϵ2\alpha ^23k`$, $`P=2ϵ3k`$, $`R=k2\alpha ^2`$, and $`\mathrm{\Gamma }=1+\delta `$. Here the width $`\alpha `$ is the only independent parameter. As the above form of the solution indicates, both nonsingular and singular solitons are possible solutions depending on the values of $`ϵ`$, and the source strength $`k`$. We have checked that periodic solutions are also allowed. Below we give a few more solutions, illustrating the differences between the localized soliton and the periodic solitary wave solutions.

Case(I): Trigonometric solution.- In the limit $`m=0`$, unlike the unperturbed NLSE, in this case, one finds rational solutions of the trigonometric type. Apart from the general solutions, interestingly, for these, one can obtain special cases, where $`A=0`$ and $`B0`$ is allowed. However, the vice-versa is forbidden. Following is an example of this type of nonsingular periodic solutions, for the repulsive case:

$$\rho (\xi )=(\frac{2k}{ϵ})\frac{\mathrm{cos}^2(\xi )}{1\frac{2}{3}\mathrm{cos}^2(\xi )}.$$
(13)
Here, $`ϵ`$ has to be negative since, $`\alpha ^2=ϵ/4`$; and is given by $`ϵ=(27gk^2/2)^{1/3}`$. This periodic solution is found to be stable, as evidenced from the numerical simulations seen in Fig.1.

Case(II): Hyperbolic solution.- Unlike the above periodic case, here one finds that, the solutions with $`B=0`$, and $`A0`$ are allowed, the vice-versa not being true. Hence, these localized solitons behave differently from the soliton trains. There are both singular and nonsingular solutions. For $`B=0`$, and $`m=1`$, we found that $`\alpha ^2=ϵ/4`$, and $`ϵ=(27gk^2/2)^{1/3}`$. This yields, the singular hyperbolic solution,

$$\rho (\xi )=(\frac{3k}{ϵ})\frac{1}{1\frac{3}{2}\mathrm{sech}^2(\xi )}.$$
(14)
The singularity here may correspond to extreme increase of the field amplitude due to self-focussing, as is known for the other nonlinear systems. We have also found nonsingular solutions of the above type.

Case(III): Pure cnoidal solutions.- In the general case, for $`0<m<1`$, from Eq. (11), one obtains the width $`\alpha `$, in terms of the other parameters. The other three equations then yield $`A`$, $`B`$, and $`D`$; in the process one encounters a cubic equation whose real roots need only to be considered. This puts restriction on the solution parameters. Interestingly, the constraints imply an $`m`$ independent condition on the parameters in the form

$`5ϵAD^2+ϵAD^3+ϵBD+5ϵBD^2+3gA^2BD`$
$`+3gADB^26kD^26kD^3+3gA^3D^2+3gB^3=0.`$
Below we list a few special cases. For $`A=0`$, $`D=1`$, and $`m=5/8`$; we found that for $`\alpha ^2=(2/7)ϵ`$, and $`ϵ=7(gk^2/18)^{1/3}`$; this corresponds to attractive case. Explicitly, the solution is given by,

$$\rho (\xi )=(\frac{14k}{3ϵ})\frac{\mathrm{cn}^2(\xi ,m)}{1+\mathrm{cn}^2(\xi ,m)}.$$
(15)
For $`A=0`$, and $`m=1/2`$; it is found that $`\alpha ^2=ϵ/2\sqrt{3}`$ and $`ϵ=(27gk^2)^{\frac{1}{3}}`$, for which

$$\rho (\xi )=(\frac{2\sqrt{3}k}{ϵ})\frac{\mathrm{cn}^2(\xi ,m)}{1+\frac{1}{\sqrt{3}}\mathrm{cn}^2(\xi ,m)}.$$
(16)
Since the localized solitons are usually robust, we have performed numerical simulations to check the stability of the solutions pertaining to Case(I), i.e., the trigonometric solution. It is worth pointing out that the numerical techniques based on the fast Fourier transform (FFT) are expensive as they require the FFT of the external source. Hence, we have used the Crank-Nicholson finite difference method press  to solve the NLSE with a source, which is quite handy, and unconditionally stable. The initial conditions chosen from the exact solution are knitted on a lattice with a grid size $`dx=0.005`$, and $`dt=5.0\times 10^6`$. The evolution of this solution, as depicted in Fig.1 indicate that it is a stable solution.

In conclusion, we have used a fractional transformation to connect the solutions of the phase-locked NLSE with the elliptic functions, in an exact manner. The solutions are necessarily of the rational type that contain solitons, solitary waves, as also singular ones. Our procedure is applicable, both for the attractive and repulsive cases. Because of their exact nature, these will provide a better starting point for the treatment of general externally driven NLSE. Considering the utility of this equation in fiber optics and other branches of physics, these solutions may find practical applications.

We acknowledge profitable discussions with Dr. E. Alam regarding the algorithm that has been implemented here, and Prof. A. Khare for many useful discussions.


text:
# Phase diagram of the frustrated Hubbard model

## Abstract

The Mott-Hubbard metal-insulator transition in the paramagnetic phase of the one-band Hubbard model has long been used to describe similar features in real materials like V2O3. Here we show that this transition is hidden inside a rather robust antiferromagnetic insulator even in the presence of comparatively strong magnetic frustration. This result raises the question of the relevance of the Mott-Hubbard metal-insulator transition for the generic phase diagram of the one-band Hubbard model.

The microscopic description of magnetism and metal-insulator transitions constitutes one of the major research activities in modern solid state theory. Especially transition metal compounds like V2O3, LaTiO3, NiS2-xSex and the cuprates show metal-insulator transitions and magnetic order depending on composition, pressure or other control parameters imada98 . The paramagnetic insulating phase observed in these materials is believed to be a so-called Mott-Hubbard insulator due to electron-electron correlations; in contrast to Slater or band insulators like SrTiO3.

The simplest model showing both magnetism and a correlation-induced metal-insulator transition (MIT) is the one-band Hubbard model hubbard

$$H=\underset{i,j,\sigma }{}t_{ij}c_{i\sigma }^{}c_{j\sigma }^{}+\frac{U}{2}\underset{i\sigma }{}n_{i\sigma }n_{i\overline{\sigma }}.$$
(1)
Considerable progress in understanding the physics of this simple but nevertheless non-trivial model has been achieved in the last decade through the development of the dynamical mean-field theory (DMFT) mv ; pradv ; rmp . In particular, the phase diagram for the unfrustrated Hubbard model is very well understood pradv ; rmp . At half-filling the physics is dominated by an antiferromagnetic insulating phase (AFI) for all $`U>0`$ with a maximum $`T_N0.15W`$ around $`UW`$, where $`W`$ is the bandwidth of the non-interacting system. For finite doping, the antiferromagnetic phase persists up to a critical doping $`\delta _c`$ jazph  and in addition shows phase separation pvd ; zitz1 . For very large values of $`U`$ the antiferromagnetic phase is replaced by a small region of Nagaoka type ferromagnetism nagoka ; oberm ; ulmke .

Under rather special assumptions muha89 , it is possible to introduce complete magnetic frustration and suppress the antiferromagnetic phase. In this case, a transition from a paramagnetic metal (PM) to a paramagnetic insulator (PI) becomes visible at half filling. At $`T=0`$ it occurs at a value of the Coulomb parameter $`U_c1.5W`$ jazph ; rmp ; bulprl . Interestingly, the transition is of first order rmp ; bulcosvol  for $`T>0`$ with a second order end point at a $`T_c0.017W`$ and $`U_c1.2W`$. Note that $`T_cT_N^{\mathrm{max}}`$.

A closer look at the phase diagram of V2O3 mcwhan  reveals a strikingly similar scenario, and indeed the DMFT results for the Hubbard model have been used as a qualitative explanation rmp ; roz95 . For a proper description of this material, however, the antiferromagnetic phase below $`T_N160`$K mcwhan  has to be taken into account. It was argued and generally accepted rmp  that the introduction of partial magnetic frustration will lead to the anticipated situation, where the MIT extends beyond the antiferromagnetic phase at low temperatures. The merging of these two transitions presents an interesting problem on its own, because it is commonly believed that the magnetic transition should be of second order. Furthermore, previous results for a system with magnetic frustration show, in addition to the desired effect of reducing $`T_N`$, an extended antiferromagnetic metallic (AFM) phase at small $`U`$, preceeding the transition to the AFI rmp ; chitra . This observation suggests an appealing possibility to link the MIT in the paramagnetic phase with a transition from an AFM to an AFI.

In this paper, we present results from a calculation using Wilson’s numerical renormalization group approach (NRG) nrg  and exact diagonalization techniques (ED) krauth  to solve the DMFT equations rmp ; bulprhew  for the Hubbard model with magnetic frustration at half filling. We show that frustration of the magnetic interactions through incorporation of suitable long-range hopping does not lead to the previously reported sequence PM$``$AFM$``$AFI with an extended region of an AFM at $`T=0`$ rmp ; chitra . Instead, we observe a first-order transition PM$``$AFI. Furthermore, the reduction of $`T_N`$ is too small to result in the qualitative phase diagram of V2O3.

The natural choice for studying the effect of magnetic frustration is the simple hypercubic lattice with nearest and next-nearest neighbor hopping. In the DMFT, the lattice structure only enters via the dispersion of the noninteracting band states, and the corresponding $`𝐤`$ sums can conveniently be transformed into energy integrals using the free density of states (DOS) rmp . A further simplification arises if one considers lattices in the limit of large coordination number. Especially for the simple hypercubic lattice the DOS then becomes a Gaussian mv ; muha89  and the integrals can be performed analytically pradv .

The investigation of magnetic properties is straightforward, too. In the case of the Néel state, the lattice is divided into A and B sublattices which results in a matrix structure of the DMFT equations rmp . An antiferromagnetic Néel order then corresponds to a finite staggered magnetization $`m_S>0`$ with $`m_A=m_S`$ and $`m_B=m_S`$. Unfortunately, the Gaussian DOS of the hypercubic lattice has no real band edges, but stretches to infinity; the resulting exponential tails therefore prevent a clear distinction between metal and insulator at $`T=0`$, as has been observed in a Hartree calculation for the hypercubic lattice with next-nearest neighbor hopping hofvol .

Georges et al. rmp  suggested an extension of the DMFT equations for the Bethe lattice which introduces magnetic frustration with the DOS having finite support; this results in an analytically tractable form of the DMFT equations even for the AB-lattice.

For the standard Bethe lattice with nearest-neighbor hopping $`t`$ and in the limit of infinite coordination number the DMFT equations on an AB-lattice acquire the form rmp

$$\begin{array}{c}G_{A\sigma }(z)=\frac{1}{z+\mu \mathrm{\Sigma }_{A\sigma }(z)\frac{t^2}{4}G_{B\sigma }(z)},\hfill \\ G_{B\sigma }(z)=\frac{1}{z+\mu \mathrm{\Sigma }_{B\sigma }(z)\frac{t^2}{4}G_{A\sigma }(z)}.\hfill \end{array}$$
(2)
Frustration is then introduced via additional terms in the denominators of (2) rmp ; chitra

$$\begin{array}{c}G_{A\sigma }(z)=\frac{1}{z+\mu \mathrm{\Sigma }_{A\sigma }(z)\frac{t_1^2}{4}G_{B\sigma }(z)\frac{t_2^2}{4}G_{A\sigma }(z)},\hfill \\ G_{B\sigma }(z)=\frac{1}{z+\mu \mathrm{\Sigma }_{B\sigma }(z)\frac{t_1^2}{4}G_{A\sigma }(z)\frac{t_2^2}{4}G_{B\sigma }(z)},\hfill \end{array}$$
(3)
and the constraint $`t^2=t_1^2+t_2^2`$. In the paramagnetic case, the equations (3) reduce to those of a standard Bethe lattice which, for example, result in the well-studied Mott-transition. Furthermore, despite the frustration introduced, the system is still particle-hole symmetric. Especially for half filling this feature reduces the numerical effort quite drastically. Note that a similar suggestion exists for the hypercubic lattice muha89 , too.

Invoking the symmetry $`G_{A\sigma }(z)=G_{B\overline{\sigma }}(z)`$ valid for the Néel state, eqs. (3) reduce to two coupled nonlinear equations which we solve iteratively. In the course of the iterations, the quantity $`\mathrm{\Sigma }_\sigma (z)`$ has to be calculated from the solution of a generalized single impurity Anderson model rmp . For that task we employ the NRG nrg , suitably extended to treat spin-polarized situations zitz1 ; hof01 .

Let us first discuss the results for the magnetization as function of $`U`$ for $`T=0`$. In the following, we fix $`t_2/t_1=1/\sqrt{3}0.58`$ comment  and use the bandwidth $`W`$ of the non-interacting system as our energy scale. The NRG results in Fig. 1 (circles)

show a completely different behavior as compared to the ED data (diamonds) from ref. chitra, . Instead of a continuous increase of the staggered magnetization $`m_S`$ for $`U>U_c0.4W`$ as suggested by both a Hartree calculation (crosses) and the data from ref. chitra, , we find a jump in $`m_S`$ at a considerably larger $`U_c0.9W`$. To clarify this discrepancy, we performed ED calculations, resulting in the triangles shown in Fig. 1. We find good agreement with our NRG results, the transition systematically approaching the NRG curve with increasing size of the system diagonalized in the ED procedure. We furthermore observed a rather strong dependence of the ED results on details of the numerical procedure, especially on the energy cutoff introduced in calculating $`G(z)`$. Decreasing this cutoff systematically shifts the ED result towards the one found in chitra . The NRG, on the other hand, is stable with respect to changes in the parameters controlling its numerical accuracy.

Another important question is the existence of an antiferromagnetic metallic solution of the DMFT equations. Figure 2 shows the NRG results for the DOS for $`T=0`$ and spin up on an A lattice site. Due to particle-hole symmetry the DOS for spin down on A sites (or spin up on B sites) can be obtained by $`\omega \omega `$. The full and dashed lines represent the AFI solution for $`UU_c`$ and the PM solution for $`UU_c`$, respectively. Clearly, the magnetic solution is insulating with a well-developed gap at the Fermi energy. Quite generally, we were not able to find a stable AFM solution at $`T=0`$.

The discontinuity in the staggered magnetization $`m_S`$ at the transition PM$``$AFI implies a first order transition and the existence of a hysteresis region.

Indeed, starting from the paramagnet at $`UU_c`$ and increasing $`U`$ results in a magnetization curve different from the one obtained by starting at $`UU_c`$ and decreasing $`U`$. This is apparent from Fig. 3 (full lines) where a region of hysteresis can be observed in the staggered magnetization (for temperature $`T=0.0155W`$). At the same time the total DOS at the Fermi energy $`A(0)=A_{}(0)+A_{}(0)`$ shows hysteresis between metallic and insulating behavior in exactly the same $`U`$ region. Note, that due to the finite temperature the DOS at the Fermi level is not exactly zero in the Néel state, but strongly reduced as compared to the metal bulcosvol .

It is of course important to verify that the hysteresis found for small $`U`$ is not some kind of artefact. This can most conveniently be shown by looking at the transition at large $`U`$. Due

to the mapping of the Hubbard model to a Heisenberg model in this regime, one should expect the transition to be of second order, with the staggered magnetization vanishing continously like $`m_S\sqrt{T_NT}`$ when approaching $`T_N`$ from below. That this is indeed the case is apparent from Fig. 4, where we show the squared staggered magnetization as function of $`T`$ for $`U/W=2`$. The transition is thus of second order with the expected mean-field exponent in this region of the phase diagram.

Collecting the results for the transitions and the hysteresis region for different temperatures leads to the phase diagram in Fig. 5.

An enlarged view of the region showing coexistence of PM and AFI is given in the inset, where the full line represents the transition PM$``$AFI with increasing $`U`$ and the dashed line the transition AFI$``$PM with decreasing $`U`$. These two lines seem to merge at a value of $`UW`$ for this particular value of $`t_2`$, with a critical temperature for this endpoint $`T_c0.02W`$. Note that, even in the presence of such a sizeable $`t_2`$, the antiferromagnetic phase still completely encompasses the paramagnetic MIT (dotted lines in the main panel of Fig. 5 bulcosvol ).

It is, of course, interesting to see how the magnetic phase evolves with increasing $`t_2`$ and in particular how its boundary crosses the paramagnetic MIT. We find that increasing $`t_2`$ does not change the form of the magnetic

phase in Fig. 5 qualitatively, but mainly shifts the lower critical $`U`$ and decreases the maximum $`T_N`$. The calculated estimates for those two quantities as function of $`t_2`$ lead to the schematic evolution of the phase diagram presented in Fig. 6a-c. Here, only the true phase boundaries are shown. A direct calculation of the free energy at finite temperatures is presently not possible with the NRG method, so we cannot calculate the actual transition line separating the paramagnetic and AF phases. The transition lines in Fig. 6a-c are therefore a guide to the eye only. For the Mott transition, the position of the actual transition line has been calculated in ref. tong, .

Figure 6a shows the qualitative phase diagram corresponding to Fig. 5 with the line of first order transitions ending in a critical point. Upon further increasing the value of $`t_2`$, the first order transition lines from both the PM$``$AFI and the Mott transition cross (Fig. 6b), thus exposing a finite region of the Mott insulator and a transition PI$``$AFI. Finally, for even higher values of $`t_2`$, the PM$``$AFI transition at $`T=0`$ approaches the Mott transition and $`T_N`$ is reduced significantly (Fig. 6c). Note that in the limiting case $`t_2=t_1`$ the AFI phase completely vanishes due to the structure of the DMFT equations (3). However, for $`t_2t_1`$ there is always a finite antiferromagnetic exchange $`J(t_1^2t_2^2)/U`$ which is sufficient to stabilize an antiferromagnetic ground state for $`U>U_c`$ of the Mott transition.

From these results we conclude that frustration as introduced via eqs. (3) is not sufficient to qualitatively reproduce the phase diagram of materials like V2O3. In particular, the Mott transition extends beyond the AFI region only for unphysically large values of $`t_2`$.

The question remains whether it is possible at all to reproduce qualitatively the scenario observed in V2O3 within a one-band model. Based on our results reported here, we rather believe that one has to take into account additional degrees of freedom, for example phonons (within a Holstein-Hubbard model) or orbital degeneracies (within a multi-band Hubbard model).

###### Acknowledgements.

We acknowledge useful conversations with M. Vojta, G. Kotliar, R. Chitra and P.J.G. van Dongen. This work was supported by the DFG through the collaborative research center SFB 484, the Leibniz Computer center and the Computer center of the Max-Planck-Gesellschaft in Garching. NT acknowledges the support by the Alexander von Humboldt foundation.


text:
# Universality in the physics of cold atoms with large scattering length11footnote 1 This work was done in collaboration with Eric Braaten.

## 1 INTRODUCTION

The Effective Field Theory (EFT) approach provides a powerful framework that exploits the separation of scales in physical systems. Only low-energy (or long-range) degrees of freedom are included explicitly, with the rest parametrized in terms of the most general local (contact) interactions. This procedure exploits the fact that a low-energy probe of momentum $`k`$ cannot resolve structures on scales smaller than $`R1/k`$. (Note that $`\mathrm{}=c=1`$ in this talk.) Using renormalization, the influence of short-distance physics on low-energy observables is captured in a small number of low-energy constants. Thus, the EFT describes universal low-energy physics independent of detailed assumptions about the short-distance dynamics. All physical observables can be described in a controlled expansion in powers of $`kl`$, where $`l`$ is the characteristic low-energy length scale of the system. The size of $`l`$ depends on the system under consideration: for a finite range potential, e.g., it is given by the range of the potential. For the systems discussed here, $`l`$ is of the order of the effective range $`r_e`$.

In this talk, I will focus on applications of EFT to few-body systems with large S-wave scattering length $`al`$. For a generic system, the scattering length is of the same order of magnitude as the low-energy length scale $`l`$. Only a very specific choice of the parameters in the underlying theory (a so-called fine tuning) will generate a large scattering length. Systems with large scattering length can be found in many areas of physics. Examples are the S-wave scattering of nucleons and of 4He atoms. For alkali atoms close to a Feshbach resonance, $`a`$ can be tuned experimentally by adjusting an external magnetic field.

## 2 THREE-BODY SYSTEM WITH LARGE SCATTERING LENGTH

In this section, we give a very brief review of the EFT for few-body systems with large scattering length $`a`$. We will focus on S-waves (For more details, see Ref. ).

For typical momenta $`k1/a`$, the EFT expansion is in powers of $`r_e/a`$ so that higher order corrections are suppressed by powers of $`r_e/a`$. The leading order corresponds to $`r_e=0`$. We consider a 2-body system of nonrelativistic bosonic atoms with large scattering length $`a`$ and mass $`m`$. At sufficiently low energies, the most general Lagrangian for S-wave interactions may be written as:

$$=\psi ^{}\left(i_t+\frac{\stackrel{}{}^2}{2m}\right)\psi \frac{C_0}{2}(\psi ^{}\psi )^2\frac{D_0}{6}(\psi ^{}\psi )^3+\mathrm{},$$
(1)
where the $`C_0`$ and $`D_0`$ are nonderivative 2- and 3-body interaction terms, respectively. The strength of the $`C_0`$ term is determined by the scattering length $`a`$, while $`D_0`$ depends on a 3-body parameter to be introduced below. The dots represent higher-order derivative terms which are suppressed at low-energies. For momenta $`k`$ of the order of the inverse scattering length $`1/a`$, the problem is nonperturbative in $`ka`$. The exact 2-body scattering amplitude can be obtained analytically by summing the so-called bubble diagrams with the $`C_0`$ interaction term. The $`D_0`$ term does not contribute to 2-body observables. After renormalization, the resulting amplitude reproduces the leading order of the well-known effective range expansion for the atom-atom scattering amplitude: $`f_{AA}(k)=(1/aik)^1,`$ where the total energy is $`E=k^2/m`$. If $`a>0`$, $`f_{AA}`$ has a pole at $`k=i/a`$ corresponding to a shallow dimer with binding energy $`B_2=1/(ma^2)`$. Higher-order derivative interactions are perturbative and give the momentum-dependent terms in the effective range expansion.

We now turn to the 3-body system. Here, it is useful to introduce an auxiliary field for the two-atom state (see Ref.  for details). At leading order, the atom-dimer scattering amplitude is given by the integral equation shown in Fig. 1.

A solid line indicates a single atom and a double line indicates an interacting two-atom state (including rescattering corrections). The integral equation contains contributions from both the 2-body and the 3-body force terms. The inhomogeneous term is given by the first two diagrams on the right-hand side: the one-atom exchange diagram and the 3-body force. The integral equation simply sums these diagrams to all orders. After projecting onto S-waves, we obtain the equation

$`𝒯(k,p;E)`$ $`=`$ $`{\displaystyle \frac{16}{3a}}M(k,p;E)+{\displaystyle \frac{4}{\pi }}{\displaystyle _0^\mathrm{\Lambda }}{\displaystyle \frac{dqq^2M(q,p;E)}{1/a+\sqrt{3q^2/4mEiϵ}}}𝒯(k,q;E),`$ (2)
for the off-shell atom-dimer scattering amplitude with the inhomogeneous term

$`M(k,p;E)`$ $`=`$ $`{\displaystyle \frac{1}{2pk}}\mathrm{ln}\left({\displaystyle \frac{p^2+pk+k^2mE}{p^2pk+k^2mE}}\right)+{\displaystyle \frac{H(\mathrm{\Lambda })}{\mathrm{\Lambda }^2}}.`$ (3)
The logarithmic term is the S-wave projected one-atom exchange, while the term proportional to $`H(\mathrm{\Lambda })`$ comes from the 3-body force. The physical atom-dimer scattering amplitude $`f_{AD}`$ is given by the solution $`𝒯`$ evaluated at the on-shell point: $`f_{AD}(k)=𝒯(k,k;E)`$ where $`E=3k^2/(4m)1/(ma^2)`$. The 3-body binding energies $`B_3`$ are given by those values of $`E`$ for which the homogeneous version of Eq. (2) has a nontrivial solution.

Note that an ultraviolet cutoff $`\mathrm{\Lambda }`$ has been introduced in (2). This cutoff is required to insure that Eq. (2) has a unique solution. All physical observables, however, must be invariant under changes of the cutoff, which determines the behavior of $`H`$ as a function of $`\mathrm{\Lambda }`$ :

$`H(\mathrm{\Lambda })={\displaystyle \frac{\mathrm{cos}[s_0\mathrm{ln}(\mathrm{\Lambda }/\mathrm{\Lambda }_{})+\mathrm{arctan}s_0]}{\mathrm{cos}[s_0\mathrm{ln}(\mathrm{\Lambda }/\mathrm{\Lambda }_{})\mathrm{arctan}s_0]}},`$ (4)
where $`s_0=1.00624`$ is a transcendental number and $`\mathrm{\Lambda }_{}`$ is a 3-body parameter introduced by dimensional transmutation. This parameter cannot be predicted by the EFT and must be taken from experiment. Note that $`H(\mathrm{\Lambda })`$ is periodic and runs on a limit cycle. When $`\mathrm{\Lambda }`$ is increased by a factor of $`\mathrm{exp}(\pi /s_0)22.7`$, $`H(\mathrm{\Lambda })`$ returns to its original value.

In summary, two parameters are required in the 3-body system at leading order in $`r_e/a`$: the scattering length $`a`$ (or the dimer binding energy $`B_2`$) and the 3-body parameter $`\mathrm{\Lambda }_{}`$ . The EFT reproduces the universal aspects of the 3-body system that were first derived by Efimov . These include the accumulation of infinitely many 3-body bound states (so-called Efimov states) at threshold as $`a\pm \mathrm{}`$. As we will demonstrate in the next section, EFT is a very efficient calculational tool to calculate those properties.

## 3 APPLICATION TO COLD ATOMS AND BEC

In this section, we discuss some applications of the EFT to systems of cold atoms and BEC’s.

First we turn to the 3-body system of 4He atoms, where $`a/r_e15`$. Both the 4He dimer and trimer were observed. The complete 3-body bound state spectrum in leading order in the EFT follows from solving Eq. (2). It can also be parametrized in terms of a universal function $`\mathrm{\Delta }`$ of the angle $`\xi =a\sqrt{mB_3}`$ and a 3-body parameter . The function $`\mathrm{\Delta }(\xi )`$ was recently calculated in EFT . Unfortunately, no 3-body observables have been measured for 4He atoms up to now, so that $`\mathrm{\Lambda }_{}`$ can not be determined from experiment. However, extensive bound state calculations with modern potentials exist. We have used recent calculations by Motovilov et al. , to obtain the dimer binding energy and the 3-body parameter $`\mathrm{\Lambda }_{}`$ from the calculated binding energy of the trimer excited state $`B_3^{(1)}`$. We have then used this input to calculate the trimer ground state energy $`B_3^{(0)}`$ for four different 4He potentials . The results are shown in Table 1.

The predicted ground state binding energy $`B_3^{(0)}`$ agrees well with the direct calculation of Motovilov et al. . This demonstrates that both the excited and ground state of the 4He trimer are Efimov states.

Universality also constrains 3-body scattering observables. For example the atom-dimer scattering length can be expressed in terms of $`a`$ and $`\mathrm{\Lambda }_{}`$ as

$$a_{AD}=a(1.462.15\mathrm{tan}[s_0\mathrm{ln}(a\mathrm{\Lambda }_{})+0.09])(1+𝒪(r_e/a)),a>0.$$
(5)
where the numerical constants were calculated in EFT. Using the value of $`\mathrm{\Lambda }_{}`$ extracted from the excited state binding energy, a good agreement with $`a_{AD}`$ from the direct calculation of Ref.  is obtained .

Universality is also manifest in universal scaling functions. In the left panel of Fig. 2, we display the scaling function relating $`B_3^{(1)}/B_2`$ to $`B_3^{(0)}/B_2`$ . The data points give various calculations using modern 4He potentials while the solid line gives the universality prediction from EFT. Different points on this line correspond to different values of $`\mathrm{\Lambda }_{}`$. The small deviations of the potential calculations from the universal curve are mainly due to effective range corrections and can be calculated at next-to-leading order in EFT. The calculation corresponding to the data point far off the universal curve can easily be identified as problematic.

Next we turn to 3-body recombination, which is the process when three atoms scatter to form a dimer and the third atom balances energy and momentum. This is one of the main loss processes for trapped atoms and condensates of atoms near a Feshbach resonance. The event rate can be parametrized as $`\nu =\alpha \rho ^3`$, where $`\rho `$ is the density of the atoms and $`\alpha `$ is the recombination constant. At threshold, the leading order EFT result for the contribution to $`\alpha `$ from recombination into the shallow dimer is

$$\alpha =67.1\frac{a^4}{m}\mathrm{sin}^2[s_0\mathrm{ln}(a\mathrm{\Lambda }_{})+0.19](1+𝒪(r_e/a)),a>0.$$
(6)
If deeply bound dimers are present, there are additional contributions from recombination into the deep bound states. If $`a<0`$, recombination can only go into the deep states.

In the right panel of Fig. 2, we show the 3-body loss rate coefficient $`K_3=3\alpha `$ from Ref.  for two different values of $`\mathrm{\Lambda }_{}`$ compared to experimental data for 3-body losses in a cold gas of 85Rb atoms near a Feshbach resonance at $`B=155`$ G . The value $`\mathrm{\Lambda }_{}=5`$/a.u. seems to be preferred by the data.

## 4 SUMMARY & OUTLOOK

The renormalization of the EFT for 3-body systems with large scattering length requires a 3-body force at leading order . The renormalization group evolution of the 3-body force is governed by a limit cycle. Via dimensional transmutation, the 3-body force introduces a dimensionful parameter, $`\mathrm{\Lambda }_{}`$, that parametrizes universal relations between different 3-body observables.

The EFT is very general and can be applied to many physical systems ranging from the 3-body system of 4He atoms , to 3-body recombination and dimer deactivation in cold atomic gases/BEC’s , to the triton and hypertriton in nuclear physics .

Future challenges include the extension of the EFT to the four-body system, the stability and phase structure of cold atomic gases/BEC’s close to a Feshbach resonance , and the possibility of coexisting condensates of atoms, dimers, and trimers .