text:
# Radiative Shock-Induced Collapse of Intergalactic Clouds

## 1 Introduction

Combined optical and radio observations of a number of extragalactic radio jets in both low-luminosity (FRI) and high-luminosity (FRII) galaxies have revealed interesting correlations between the jets and apparent regions of active star formation. One of the first objects demonstrated to show such a correlation was the nearest radio galaxy, Centaurus A (e.g. Blanco et al., 1975). Other examples have been found as the sensitivity and spatial resolution of radio and optical telescopes has improved. In the case of Centaurus A, recent observations with the Hubble Space Telescope, when compared to radio images obtained with the Very Large Array, have confirmed that there are about half a dozen young ($`<15`$ Myr) OB associations near filaments of ionized gas located between the radio jet and a large H I cloud (Mould et al., 2000; Rejkuba et al., 2002). Another nearby example is “Minkowski’s Object,” a peculiar small starburst system at the end of a radio jet emanating from the elliptical galaxy NGC 541, located near the center of the cluster of galaxies Abell 194 (van Breugel et al., 1985). Star forming regions associated with radio sources have also been found in cooling flow clusters (McNamara, 2002).

Correlations between radio and optical emissions have also been observed in the so-called “alignment effect” in distant ($`z>0.6`$) radio galaxies (Chambers, Miley, & van Breugel, 1987; McCarthy et al., 1987). The best studied example here is 4C41.17 at $`z=3.8`$, where deep spectroscopic observations have shown that the bright, spatially extended, rest-frame UV continuum emission aligned with the radio axis of this galaxy is unpolarized and shows P Cygni-like features similar to those seen in star-forming galaxies (Dey et al., 1997).

Collectively, these observations are most compellingly explained by models in which shocks generated by the radio jet propagate through an inhomogeneous medium and trigger gravitational collapse in relatively overdense regions (Begelman & Cioffi, 1989; De Young, 1989; Rees, 1989). A detailed analysis of the jet-induced star formation in 4C41.17 has been presented by Bicknell et al. (2000). In that object, Hubble Space Telescope images showed a bimodal optical continuum structure parallel to the radio jet (van Breugel et al., 1999), strongly supporting the idea that the star formation was triggered by sideways shocks.

Although the analytic arguments in favor of this model are compelling, the complex nature of this problem, including nonlinear effects from higher-order coupling between hydrodynamics, cooling, and gravity, suggest that numerical studies may provide a more complete and detailed picture. Furthermore, numerical studies allow us to probe environments that may be observationally out of reach. For example, such studies may be able to parameterize the conditions of shock-induced star formation in the early universe, where direct observations of star-forming regions are not possible. This could be important to understanding the role of jets in the feedback of active galactic nuclei (AGN) on their environment. Numerical simulations also give us insight into the role of shock-induced star formation in other environments such as supernova shocks and cloud-cloud collisions.

The interaction of a strong shock with a single, non-radiative cloud has been the subject of many numerical studies. A thorough analysis of this problem and a review of relevant literature is provided by Klein, McKee, & Colella (1994). A more recent study, considering the interaction of a strong shock with a system of clouds, is presented by Poludnenko, Frank, & Blackman (2002). For a non-radiative cloud the passing shock ultimately destroys the cloud within a few dynamical timescales. Destruction results primarily from hydrodynamic instabilities at the interface between the cloud and the post-shock background gas.

Less well studied is the case of a strong shock interacting with a radiative cloud. Mellema, Kurk, & Röttgering (2002) showed that the effects of a shock on such a cloud are very different from the non-radiative case. Instead of re-expanding and quickly diffusing into the background gas, the compressed cloud instead breaks up into numerous dense, cold fragments. These fragments survive for many dynamical timescales and are presumably the precursors to star formation.

In this work, we extend previous numerical studies by considering two primary cooling models. The first model uses an equilibrium cooling curve with a tunable metallicity, set to solar in this work. As such, this model provides a reasonable approximation for cooling in metal-enriched clouds around low-redshift (FRI) radio galaxies. The second model solves the full non-equilibrium chemistry for a primordial gas mixture. This model is very well suited to address jet-induced star formation in pristine clouds around high-redshift (FRII) galaxies. With both models we explore a substantial subset of the parameter space over which shock-induced star formation may apply. We proceed in § 2 by discussing the relevant physical processes, the timescales over which they act, and the limitations of the present study. The numerical methods and models employed in the current work are discussed in § 3, and the results of the models are presented in § 4. Further implications for star formation are discussed in § 5. Our conclusions are briefly recapped in § 6.

## 2 Shock-Cloud Interactions

The basic idea of shock-induced star-formation is relatively simple. A strong shock passes through a clumpy medium, triggering many smaller-scale compressive shocks in overdense clumps. These shocks increase the density inside the clumps and make it possible for them to radiate more efficiently. If the radiative efficiency of the gas has a sufficiently shallow dependence upon the temperature, then radiative emissions are able to cool the gas rapidly, in a runaway process, producing even higher densities as the cooling gas attempts to re-attain pressure balance with the surrounding medium (Field, 1965; Murray & Lin, 1989). Both the reduction in temperature and the increase in density act to reduce the Jeans mass, above which gravitational forces become important. Any clump that was initially close to this instability limit will be pushed over the edge and forced into gravitational collapse. Thus, the passing of a shock through a clumpy medium may trigger a burst of star formation. Several processes, acting on different timescales as discussed below, govern whether or not cooling and star formation are able to proceed.

### 2.1 Jet-Cloud Collision

In this work we are primarily interested in the interactions occurring over a relatively small region of a much larger intergalactic environment. However, since this work is being presented in the context of jet-induced star formation we wish to first motivate our treatment of the problem as either a single or a few overdense clumps being overrun by a planar shock. The presumption that the intergalactic medium is clumpy and multi-phase with a range of densities and temperatures is consistent with observational and theoretical studies of such objects as high-pressure ($`10^5`$ K cm-3) cooling flow clusters (Ferland, Fabian, & Johnstone, 2002) and the giant emission-line halos found around distant radio galaxies (Reuland et al., 2003). The planar-shock assumption is also consistent, although the specific relation between the shock and the jet may vary somewhat from system to system.

For luminous FRII galaxies, with their powerful relativistic jets, it is unlikely that star formation will occur within the stream of the jet itself. More likely, star formation will proceed within a cocoon of shock gas around the jet (Begelman & Cioffi, 1989; De Young, 1989; Rees, 1989). Because of the very large sizes of these cocoons, the triggering shocks will appear to be nearly planar over the length scales considered in this study.

For low-luminosity FRI galaxies, the jets have a lower Mach number and appear turbulent; no hotspots or cocoons form (Bicknell, 1984). Nevertheless, such jets may still trigger potential star-forming shocks through collisions with nearby clouds, as illustrated in the following numerical simulation. This simulation is performed using the Cosmos code described in § 3 below. Our parameters are chosen to closely resemble the conditions of Minkowski’s Object, the small starburst system at the end of the radio jet emanating from the radio galaxy NGC 541. Specifically, we set up an elliptical cloud with a major axis of 10 kpc, a minor axis of 5 kpc, a temperature $`T=10^6`$ K, and density $`n_{cl}=0.1`$ cm-3. This corresponds to an initial cloud mass of $`10^9M_{}`$.

The radio source associated with NGC 541 is of relatively low luminosity and exhibits a FRI-type morphology (Fanaroff & Riley, 1974). We, therefore, used the detailed radio and X-ray study of the proto-typical FRI-type radio galaxy 3C31 by Laing & Bridle (2002) to estimate plausible jet parameters near Minkowski’s Object, at $`15`$ kpc from the NGC 541 AGN. For our numerical simulations, we assume that the long axis of the cloud is aligned with a jet of high velocity ($`9\times 10^4`$ km s-1), low density ($`10^4`$ cm-3), hot ($`10^9`$ K) gas flowing onto our computational grid. The diameter of the jet nozzle is equal to half the diameter of the cloud along the minor axis. The assumption of an elliptical cloud is not essential to this discussion, although it can be easily motivated by tidal stretching from the nearby galaxies. It is also not essential that the jet be aligned with the long axis of the cloud, although this is convenient for numerical purposes.

As shown in Figure 1 the jet-cloud collision triggers a nearly planar shock down the long axis of the cloud. As the bow shock from the jet wraps around the cloud, it also triggers shocks along the sides of the cloud. A similar shock structure may explain the filamentary nature of the star-forming region in Minkowski’s Object. Note the similarities between the next-to-last frame of the numerical simulations and the observations as illustrated in Figure 2.

We stress that this simulation is only intended to motivate our subsequent treatment of the problem at much smaller physical scales. This simulation is unable to resolve the many moderately dense ($`10`$ cm-3), warm ($`10^4`$ K) sub-kiloparsec clouds presumed to be interspersed within the larger cloud simulated here. Further, the jet is not fully evolved when it hits the cloud, although this may not be critical since even a fully-evolved jet would trigger shocks in any clouds that crossed its path.

### 2.2 Timescales

The initial conditions for all of our subsequent calculations consist of a spherical cloud of radius $`R_{cl}`$ and density $`\rho _{cl}=m_Hn_{cl}`$ in initial pressure equilibrium with a background gas at temperature $`T_{b,i}`$ and density $`\rho _{b,i}=m_Hn_{b,i}=\rho _{cl}/\chi <\rho _{cl}`$, where $`\chi >1`$ is the density ratio between the cloud and background gases. Given this simplified geometry and parameterization, a number of useful dynamical timescales can be estimated analytically in order to weigh the relative importance of various physical effects including: shock crossing, shock compression, instability growth, cooling, freefall, magnetic acceleration, and thermal conduction.

The first timescale of interest is the shock-passage time. A planar shock of velocity $`v_{sh,b}`$ measured in the background gas, passes over the cloud in a time

$$t_{sp}=\frac{2R_{cl}}{v_{sh,b}}.$$
(1)
The gas in the cloud will react to the shock on roughly a sound-crossing time

$$t_{sc}=\frac{2R_{cl}}{c_{s,cl}},$$
(2)
where $`c_{s,cl}`$ is the initial isothermal sound speed of the gas in the cloud. Provided $`c_{s,cl}v_{sh,b}`$, the cloud will rapidly find itself in an overpressure region and a nearly spherical shock will be driven into the cloud. The velocity of this shock can be estimated by assuming pressure equilibrium between the post-shock cloud material and the shocked background region. If the shock in the background is strong ($`=v_{sh,b}/c_{s,b}1`$), then the post-shock pressure is approximately equal to $`\rho _{b,i}v_{sh,b}^2`$. Similarly, the post-shock pressure in the cloud is of order $`\rho _{cl}v_{sh,cl}^2`$. Equating these, we estimate the velocity of the shock in the cloud as

$$v_{sh,cl}\left(\frac{\rho _{b,i}}{\rho _{cl}}\right)^{1/2}v_{sh,b}=\frac{v_{sh,b}}{\chi ^{1/2}}.$$
(3)
\[More rigorous estimates are presented in Sgro (1975) and Klein, McKee, & Colella (1994).\] This shock will compress the cloud on a timescale of approximately

$$t_{cc}=\frac{R_{cl}}{v_{sh,cl}}\chi ^{1/2}\frac{R_{cl}}{v_{sh,b}}.$$
(4)
This is an important dynamical timescale for the cloud.

After the shock has passed over it, the cloud becomes susceptible to both Rayleigh-Taylor and Kelvin-Helmholtz instabilities. The growth of the Rayleigh-Taylor instability is due to the acceleration of the cloud by the post-shock background gas. The acceleration timescale is the amount of time it takes to accelerate the cloud gas to the velocity of the post-shock background. If we consider the momentum transfer from a column of post-shock background gas with velocity $`v_{b,ps}`$ and cross-sectional area $`\pi R_{cl}^2`$ to a spherical cloud, then the acceleration timescale is given roughly by $`t_{acc}\chi R_{cl}/(3v_{b,ps})`$, where we have assumed a mean post-shock force acting on the cloud of $`\pi R_{cl}^2\rho _{b,ps}v_{b,ps}^2`$ and used the shock-jump condition $`\rho _{b,ps}=4\rho _{b,i}`$ with an adiabatic index of $`\mathrm{\Gamma }=5/3`$. The shock-jump conditions also give us $`v_{b,ps}=3v_{sh,b}/4`$. Thus, the acceleration timescale can be rewritten as

$$t_{acc}\chi \frac{R_{cl}}{v_{sh,b}}\chi ^{1/2}t_{cc},$$
(5)
where we have dropped a numerical factor of order unity. This gives an acceleration of $`gv_{sh,b}/t_{acc}R_{cl}/t_{cc}^2`$, corresponding to a Rayleigh-Taylor growth time of $`t_{RT}^1(gk)^{1/2}`$ for a perturbation of wavenumber $`k`$ (Chandrasekhar, 1961). Written in terms of the cloud-compression timescale, $`t_{cc}`$, the Rayleigh-Taylor growth time is

$$t_{RT}\frac{t_{cc}}{(kR_{cl})^{1/2}}.$$
(6)
Although the shortest wavelengths ($`\lambda R_{cl}`$) have the fastest growth rates, they also saturate rapidly when their amplitudes $`A\lambda R_{cl}`$. Wavelengths corresponding to $`kR_{cl}1`$ are therefore the most disruptive. Hence, Rayleigh-Taylor instabilities will break up the cloud over a timescale comparable to the cloud compression.

For $`\chi 1`$ the timescale for growth of the Kelvin-Helmholtz instability is $`t_{KH}^1=kv_{rel}/\chi ^{1/2}`$ (Chandrasekhar, 1961), where $`v_{rel}`$ is the relative velocity between the post-shock background and the cloud. Since the cloud accelerates rather slowly for $`\chi 1`$, the relative velocity is approximately equal to the post-shock velocity of the background gas $`v_{b,ps}=3v_{sh,b}/4`$. Again, in terms of the cloud-compression timescale $`t_{cc}`$ , we have

$$t_{KH}\frac{t_{cc}}{(kR_{cl})}.$$
(7)
The Kelvin-Helmholtz instability will thus also act on a cloud-compression timescale.

The timescales above have been derived assuming adiabatic evolution of the gas and would be modified somewhat if radiative cooling were taken into account. However, our primary purpose is to roughly characterize some general timescales for these clouds, for which the values derived above are more than adequate.

If the cloud is not able to cool radiatively, then the above results suggest that the cloud will be destroyed over a time comparable to the cloud-compression timescale. This has been shown to indeed be the case in numerous hydrodynamic simulations (e.g. Klein, McKee, & Colella, 1994; Poludnenko, Frank, & Blackman, 2002). Destruction is enhanced by the fact that, after the compressive shock passes through itself and exits at the opposite surface of the cloud, it triggers a rarefaction wave traveling back into the cloud. This causes the cloud to re-expand and increases the pressure contrast between the compressed cloud and the post-shock external gas, effectively accelerating the growth of destructive instabilities.

However, if the cooling timescale is short compared to the cloud-compression timescale ($`t_{cool}t_{cc}`$), a very different outcome is reached, as was shown in Mellema, Kurk, & Röttgering (2002). Cooling instabilities in the post-shock cloud material cause it to collapse into a very thin, cold dense shell behind the shock. This greatly enhances the density contrast between the post-shock cloud and background gases, which slows the growth of destructive instabilities. We can estimate the cooling timescale as

$$t_{cool}=\frac{1.5n_{cl,ps}k_BT_{cl,ps}}{(n_{cl,ps})^2\mathrm{\Lambda }(T)},$$
(8)
where the numerator is the energy density of the gas and the denominator is the volume emissivity. To get a general estimate of the cooling timescale, we use the approximate cooling rate, $`\mathrm{\Lambda }=1.33\times 10^{19}T^{1/2}`$ erg cm3 s-1 from Kahn (1976), appropriate for gas temperatures ranging from $`5\times 10^45\times 10^7`$ K. Applying the shock-jump conditions, the cooling timescale can be rewritten as

$$t_{cool}=C\frac{v_{sh,b}^3}{\chi ^{3/2}\rho _{cl}},$$
(9)
where $`C=7.0\times 10^{35}`$ g cm-6 s4. In this study we will predominantly explore the cooling-dominated regime $`t_{cool}<0.1t_{cc}`$. (A version of this condition is given in Mellema, Kurk, & Röttgering (2002), Eq. 1.) For our typical data with cloud radius $`R_{cl}=100`$ pc and density ratio $`\chi =10^3`$, we find that cooling will generally govern evolution for moderate cloud densities and shock velocities. However, sufficiently high shock velocities can suppress the effect of cooling over the cloud destruction time (e.g., for shock velocities $`v_{sh,b}>10^4`$ km s-1 with $`n_{cl}1`$ cm-3). Cooling is also negligible at sufficiently low densities (e.g., clouds with $`n_{cl}<10^4`$ cm-3 do not cool for $`v_{sh,b}10^3`$ km s-1).

### 2.3 Limitations

Our models are subject to some limitations. Due to resolution requirements and computational limitations, the models have been studied in two-dimensional, Cartesian geometry, and so the clouds represent slices through infinite cylinders. As can be seen in the models presented below, the initial compressions are highly symmetric, and so the additional convergence expected in three-dimensional models might well lead to stronger compressions and enhanced cooling, especially for the primordial clouds which rely upon the formation of H2. Three-dimensional simulations would also provide an additional degree of freedom for fragmentation through dynamical instabilities. Such models might, therefore, be expected to lead to the formation of greater numbers of fragments than seen below.

In all but one of our models, self-gravity is neglected, due to the geometry; self-gravity in cylindrical symmetry is noticeably different than in spherical symmetry. The relevant timescale for self-gravity is the free-fall timescale of the gas, $`t_{ff}=(3\pi /32G\rho )^{1/2}`$. The initial overdense clumps considered here are specifically chosen not to be self-gravitating, and the corresponding freefall time is long compared to other relevant timescales. However, as the clouds compress, the local free-fall timescale becomes significantly shorter. For most of the runs, self-gravity is just becoming important around the time we stop the simulations. In order to explore the effects of self-gravity, we include one run which solves for the self-gravity of the cloud. The results of that model indicate that the inclusion of self-gravity would likely enhance compressions and cooling in all of the models.

Examining three dimensional models with self-gravity, at comparable resolution to that used in the models presented here, will require the implementation of adaptive mesh refinement to concentrate resolution only over those regions which form dense fragments. That capability is currently being added to the code described below, and results will be presented in future work.

If the background medium is magnetized, then magnetic field lines can become trapped in deformations on the surface of the cloud. As these field lines are stretched, the magnetic pressure along the leading edge of the cloud can increase enough to accelerate the disruption of the cloud through the Rayleigh-Taylor instability (Gregori et al., 1999). We can estimate the importance of this effect by accounting for the increased acceleration of the cloud due to magnetic pressure (Gregori et al., 1999)

$$g=\frac{\rho _{b,i}v_{sh,b}^2+P_B}{R_{cl}\rho _{cl}},$$
(10)
where $`P_B=B^2/8\pi P_{b,ps}/\beta `$. The Rayleigh-Taylor growth time $`t_{RT}^1(gk)^{1/2}`$ now becomes

$$t_{RT}\frac{t_{cc}}{\{[1+3/(4\beta )]kR_{cl}\}^{1/2}}.$$
(11)
From this estimate we see that, even for strong magnetic fields ($`\beta 10`$), the Rayleigh-Taylor growth time remains comparable to the cloud-compression timescale and hence longer than the cooling timescales for most of the runs considered.

Tangled magnetic fields within the clouds would also act to resist compression, potentially reducing cooling, and enhancing cloud destruction by shocks. We are currently working to add MHD capabilities to our numerical code, and shall examine the effects of realistic magnetic field strengths in future work.

If thermal conduction is important then the cloud will evaporate into the hot background medium. The rate at which the cloud evaporates can be written (Klein, McKee, & Colella, 1994)

$$\dot{M}_{cl}=4\pi R_{cl}^2\rho c_sF(\sigma _0),$$
(12)
where $`\sigma _0`$ is the saturation parameter and $`F(\sigma _0)`$ is a dimensionless quantity of order unity. If we define an ablation timescale $`t_{ab}M_{cl}/\dot{M}_{cl}`$, then in the post-shock environment

$$t_{ab}\frac{\chi ^{1/2}t_{cc}}{7F(\sigma _0)}.$$
(13)
We see that for $`\chi 10^2`$, the ablation timescale will be comparable to or longer than the compression timescale. Thus, for the parameters considered in this study, thermal conduction can safely be ignored.

## 3 Numerical Methods

The numerical calculations discussed below have been computed using Cosmos, a massively parallel, multidimensional, radiation-chemo-hydrodynamics code designed for both Newtonian and relativistic flows developed at Lawrence Livermore National Laboratory. The relativistic capabilities and tests of Cosmos are discussed in Anninos & Fragile (2003). Tests of the Newtonian hydrodynamics options and of the microphysics relevant to the current work are presented in Anninos, Fragile & Murray (2003), and so we shall not discuss those in detail here. The calculations are carried out on a fixed, two-dimensional Cartesian ($`x`$,$`y`$) grid, implying that the simulated clouds are cylindrical rather than spherical. This limitation is currently necessary in order to maintain sufficient resolution to follow the fragmentation of the clouds.

### 3.1 Cooling Models

All of the remaining results presented here were performed at a fixed spatial resolution of $`0.5\times 0.5`$ pc per zone. There are 200 zones across the initial radius of the cloud, consistent with the resolution requirements suggested by Klein, McKee, & Colella (1994). We note, however, that in the presence of cooling, which leads to even higher densities, the resolution requirements become even more stringent. We find that we are only able to reliably follow the fragmentation of the cloud for slightly longer than the cloud-compression timescale, $`t_{cc}`$. Beyond that point, further compression of the cloud is prevented by numerical resolution rather than any physical mechanism. We have allowed a few of our runs to evolve for longer times (on the order of a few $`t_{cc}`$) to confirm that the cold, dense cloud fragments which form are relatively long lived, consistent with the findings of Mellema, Kurk, & Röttgering (2002). Further study of the evolution of these dense fragments will require finer resolution or an adaptive mesh.

#### 3.1.1 Equilibrium Cooling Curve - Low Redshift Systems

Two radiative cooling and heating models are considered in this study. In the first model, appropriate for enriched clouds around low-redshift galaxies, local cooling is given by the following cooling function:

$$\mathrm{\Lambda }(T,n)=\left[\underset{i}{}\dot{e}_i(T)(f_In)^2+\dot{e}_M(T)f_Mn^2+Jn\right]\times \{\begin{array}{cc}\mathrm{exp}[(TT_{min})/\delta T]\hfill & \text{if}TT_{min},\hfill \\ 1\hfill & \text{otherwise},\hfill \end{array}$$
(14)
based in part on an equilibrium (hydrogen recombination and collisional excitation) cooling curve. Here $`\dot{e}_i`$ is the cooling rate from hydrogen and helium lines, $`\dot{e}_M`$ is the temperature-dependent cooling rate from metals (including carbon, oxygen, neon, and iron lines, assuming solar metallicity), $`f_M`$ is a weighting or tracer function for metal cooling taken to be unity, $`J`$ is a generic background heating term, $`f_I`$ is an estimate of the ionization fraction, defined as $`\text{min}(1,\text{max}(0,(T_{eV}T_c)/3))`$ with $`T_c=1eV`$ to match roughly the expected upper and lower bounds in a mostly hydrogen gas equilibrium model, $`T=(\mathrm{\Gamma }1)e\mu /(k_Bn(1+f_I))`$ is the gas temperature in Kelvin, $`n`$ is the number density of the gas, $`e`$ is the internal energy density of the gas, $`\mu `$ is the mean molecular weight, assumed to be unity, and $`k_B`$ is Boltzmann’s constant. The exponential, with width $`\delta T=1`$ K, is introduced to suppress cooling below $`T_{min}=10`$ K. The cooling rate for metals, $`\dot{e}_M`$, is extended to low temperatures ($`10`$ K) using the curves of Dalgarno & McCray (1972) in the low ionization limit. The background heating ($`J`$) is set to initially balance the cooling inside the cloud and remains fixed throughout the evolution, consistent with the treatment of Mellema, Kurk, & Röttgering (2002).

#### 3.1.2 Non-equilibrium Primordial Chemistry - High Redshift Systems

The second cooling model applies when the chemistry is solved dynamically with the full non-equilibrium equations as described in Anninos, Fragile & Murray (2003). Currently, nine atomic and molecular species are included in the chemistry model: H I, H II, He I, He II, He III, $`e^{}`$, H-, H2, H$`{}_{}{}^{+}{}_{2}{}^{}`$. A total of 27 gas-phase chemical reactions are included in the full network. As such, this model is useful for considering pristine clouds around high redshift galaxies. The various ionization states and concentration densities $`n_i`$ of each species are calculated from the time-dependent chemistry equations using a sequential backwards differencing scheme (Anninos et al., 1997) and used explicitly in the cooling function as

$$\mathrm{\Lambda }(T,n^{[m]})=\underset{i=1}{\overset{N_s}{}}\underset{j=1}{\overset{N_s}{}}\dot{e}_{ij}(T)n^{[i]}n^{[j]}+\underset{i=1}{\overset{N_s}{}}J_in^{[i]}+\dot{e}_M(T)f_Mn^2,$$
(15)
where $`\dot{e}_{ij}(T)`$ are the cooling rates from 2-body interactions between species $`i`$ and $`j`$, and $`J_i`$ represents frequency-integrated photoionization and dissociation heating. The equation of state for temperature in this case is given by $`T=(\mathrm{\Gamma }1)e/(k_B_in_i)`$. We account for a total of seven different cooling and heating mechanisms: collisional-excitation, collisional-ionization, recombination, bremsstrahlung, metal-line cooling (dominantly carbon, oxygen, neon, and iron), molecular-hydrogen cooling, and photoionization heating. The models with non-equilibrium chemistry also include photoionization of H and photodissociation of H2 by free-streaming radiation fields. For most models, the photoionization rate is taken to be the appropriate rate expected from cosmic UV background radiation at low redshift (Bechtold et al., 1987), while the photodissociation rate is taken to be the value applicable to the local interstellar medium (Spaans & Neufeld, 1997). We also examine excursions by an order of magnitude from those values.

The essential difference between this cooling model from the equilibrium model of section 3.1.1 is that the atomic and molecular reactions are properly taken into account in order to resolve temporal phase differences between the cooling and recombination times. This allows us to more accurately predict the concentration of residual free electrons and ions from the slower recombination processes, particularly at and below the hydrogen-line cooling curves which dominate cooling down to about $`10^4`$K. These residual electrons are captured predominately by neutral hydrogen to form H-, which subsequently produces hydrogen molecules through collisional interactions with other hydrogen atoms. Cooling below the hydrogen Lyman-$`\alpha `$ line edge is achieved through the excitation of the vibrational/rotational modes of hydrogen molecules, provided H2 forms in sufficient abundance. For molecular line excitations, we use the cooling function of Lepp & Shull (1983), from which we can expect additional cooling down to about a couple hundred degrees Kelvin for our typical cloud parameters and limited grid resolution. Presumably the inclusion of dust grain physics or deuterium and higher order chemistry would contribute to further cooling. We will investigate these effects in future papers, but in this present work, we do not expect to achieve the same level of cooling in the chemistry model as in the equilibrium model, since the latter includes low temperature ($`<100`$K) cooling from metals.

In the non-equilibrium chemistry model, we also account for the self-shielding of the cloud against photoionizing and photodissociating background fields. We do this by calculating an optical depth

$$\tau _i^j(\mathrm{})=_\mathrm{}_0^{\mathrm{}}\sigma _in_i𝑑\mathrm{}^{}$$
(16)
for the two species of interest: H I and H2. The limits of integration run from the edge of the grid to points within the cloud, along the four cardinal directions (see below). The cross sections used are $`\sigma _{\mathrm{HI}}=6.3\times 10^{18}`$ cm2 (Osterbrock, 1989) and $`\sigma _{\mathrm{H}_2}=5.2\times 10^{18}`$ cm2 (Hollenbach, Werner, & Salpeter, 1971). These are the peak values of the respective cross sections. The photoionization cross section of H drops off rapidly above the ionization threshold, and so photoionization is dominated by photons near 13.6 eV, making that the most appropriate value for our approximate treatment. Photodissociation of H2 occurs for a narrow range of energies below 13.6 eV, and so the cross section can be considered to be constant. For high column densities, the wings of the absorption become important, for which the cross section would be less than used above. Our models do not achieve extremely high cross sections, and so we use the larger value of $`\sigma _{\mathrm{H}_2}`$. The optical depth calculation is initiated from each of the outer boundary cells of the computational grid and integrated along the Cartesian axes perpendicular to the boundary cell faces. For any interior zone located at position $`x`$, the optical depth is taken to be the minimum of all optical depth integrations which intersect that cell position from every exterior boundary. Hence we take

$$\tau _i(x)=\text{min}[\tau _i^x(\mathrm{}),\tau _i^x(\mathrm{}),\tau _i^y(\mathrm{}),\tau _i^y(\mathrm{})],$$
(17)
where the superscript indicates the line-of-sight direction of integration intersecting at a cell center $`x`$. Any symmetry boundaries, such as the lower $`y`$-boundary in this work, are not included in this calculation, since the optical depth in this direction would necessarily be higher. The appropriate reaction rates and heating coefficient are modified to account for the integrated absorption as

$`k20`$ $``$ $`k20\mathrm{exp}(\tau _{\mathrm{HI}}),`$
$`k27`$ $``$ $`k27\mathrm{exp}(\tau _{\mathrm{H}_2}),`$
$`J_{\mathrm{HI}}`$ $``$ $`J_{\mathrm{HI}}\mathrm{exp}(\tau _{\mathrm{HI}}),`$ (18)
following the notation of Anninos, Fragile & Murray (2003).

### 3.2 Parameter Space

For all of our calculations we have considered $`R_{cl}=100`$ pc, $`T_{cl}=10^4`$ K, and $`\chi =10^3`$ to be fixed. Assuming that the background is initially in pressure equilibrium with the cloud, we are left with only a two-dimensional parameter space to explore, $`n_{cl}`$ vs. $`v_{sh,b}`$, for the equilibrium cooling model. For the non-equilibrium chemistry model we will also consider additional parameters related specifically to that cooling model. For this combination of fixed parameters, the sound-crossing time for the cloud in every case is $`t_{sc}=17`$ Myr.

Before conducting any simulations, we first wish to narrow the parameter space to be studied. Since we are primarily interested in shock-induced star formation, we want to concentrate on clouds that are able to cool efficiently during their compression. We can define an approximate cooling-dominated regime from the timescales for cooling and cloud-compression: $`t_{cool}<0.1t_{cc}`$. Rewriting this in terms of the variables used in our numerical study, we find that cooling dominates provided that

$$\left(\frac{R_{cl}}{100\mathrm{p}\mathrm{c}}\right)>1.2\times 10^4\left(\frac{\chi }{10^3}\right)^2\left(\frac{v_{sh,b}}{10^3\mathrm{km}/\mathrm{s}}\right)^4\left(\frac{n_{cl}}{\mathrm{cm}^3}\right)^1.$$
(19)
(A similar condition is given in Mellema, Kurk, & Röttgering (2002), Eq. 1.) In the context of shocks driven by jets from AGN, this illustrates that cooling can dominate over much of the parameter space of interest. The dashed line in Figure 3 separates the cooling-dominated regions from the non-cooling for a fixed cloud radius of $`R_{cl}=100`$ pc and density ratio $`\chi =10^3`$.

We arrive at another constraint by considering the free-fall timescale of the radiatively cooled gas, $`\tau _{ff}`$, which represents the shortest possible timescale on which stars can form within the cold gas. It also represents a typical spread over which star formation shall occur. As the first massive stars form, they provide an extra heating source in the cloud that will likely shut off any further star formation (Lin & Murray, 2000; Dong, Lin, & Murray, 2003). If we take 1 Myr as a typical star-formation timescale, then any gas with a longer free-fall timescale is unlikely to form stars efficiently. Thus we define $`t_{ff}1`$ Myr as our criterion for efficient star-formation. We can estimate the free-fall timescale $`t_{ff}=(3\pi /32G\rho )^{1/2}`$ by assuming that the cold cloud gas will eventually reach pressure equilibrium with the post-shock background. With this assumption, the final density of the cloud material is

$$\rho _{cl,f}\frac{\mu m_HP_{b,ps}}{k_BT_f}.$$
(20)
Using the shock-jump condition $`P_{b,ps}\rho _{b,i}v_{sh,b}^2=(\rho _{cl}v_{sh,b}^2)/\chi `$, we can write our star-formation criterion as

$$\left(\frac{t_{ff}}{\mathrm{Myr}}\right)1.5\left(\frac{\chi }{10^3}\right)^{1/2}\left(\frac{n_{cl}}{\mathrm{cm}^3}\right)^{1/2}\left(\frac{v_{sh,b}}{10^3\mathrm{km}/\mathrm{s}}\right)^1\left(\frac{T_f}{100\mathrm{K}}\right)^{1/2}.$$
(21)
The dot-dashed line in Figure 3 shows this star-formation cut-off for specific values of $`R_{cl}=100`$ pc, $`\chi =10^3`$, $`T_f=100`$ K, and $`t_{ff}=1`$ Myr.

Two additional constraints come from the requirements that the cloud not be Jean’s unstable initially and that the shock velocity necessarily exceed the sound speed in the background gas. The first requirement gives us a constraint relation between the radius, temperature, and density of the cloud

$$\left(\frac{n_{cl}}{\mathrm{cm}^3}\right)<91\left(\frac{R_{cl}}{100\mathrm{p}\mathrm{c}}\right)^2\left(\frac{T_{cl}}{10^4\mathrm{K}}\right).$$
(22)
The horizontal solid line in Figure 3 shows this stability limit for $`R_{cl}=100`$ pc and $`T_{cl}=10^4`$ K. The second requirement can be written as a relation between the velocity of the shock and the temperature of the cloud

$$\left(\frac{v_{sh,b}}{10^3\mathrm{km}/\mathrm{s}}\right)>3.7\times 10^2\left(\frac{\chi }{10^3}\right)^{1/2}\left(\frac{T_{cl}}{10^4\mathrm{K}}\right)^{1/2}.$$
(23)
This constraint is easily met by all of the parameters considered in this study and is not shown in Figure 3 (it would be a vertical line just off the left edge of the figure).

Thus we can divide our parameter space into five regions as shown in Figure 3. Region I (gravitationally unstable) is not of interest in this study. In regions IV and V cooling is expected to be negligible, so these are not of much interest either. Although cooling will be important in region III, the star-formation timescale is too long for the process to be very efficient. This leaves us with region II as the sector of primary interest in this study of shock-induced star formation. However, we also examine models in region V in order to confirm our above estimates of the location of the boundary between regions. The particular parameter pairs chosen for numerical study (marked with crosses in Figure 3) are discussed below.

## 4 Results

Table 1 summarizes the runs using the equilibrium cooling-curve model with solar metallicities. We choose ten combinations of $`n_{cl}`$ and $`v_{sh,b}`$ in order to explore the parameter space of interest indicated in Figure 3. The Mach number, $``$, corresponding to each value of $`v_{sh,b}`$ is also listed. Finally, Table 1 also includes the shock-passage ($`t_{sp}`$), cloud-compression ($`t_{cc}`$), and cooling ($`t_{cool}`$) timescales for each run. Notice that the shock-passage and cloud-compression timescales vary by less than an order of magnitude, whereas the cooling timescale varies by almost 5 orders of magnitude among the models considered. This again arises from the very strong dependence of $`t_{cool}`$ on $`v_{sh,b}`$ and $`n_{cl}`$.

Table 2 summarizes the runs using the non-equilibrium primordial chemistry model. We choose a subset of six parameter combinations from those considered with the equilibrium cooling-curve model. Table 2 also includes the shock-passage ($`t_{sp}`$) and cloud-compression ($`t_{cc}`$) timescales for each run.

Figure 4 illustrates a sequence of density contour plots taken from the simulation of model E3. The parameters for this model were intentionally chosen to match those of run A in Mellema, Kurk, & Röttgering (2002). As expected, our results match theirs quite well. The strong cooling behind the compression shock in the cloud causes a very dense shell to form. The density of this shell continues to increase as it follows the compression shock into the center of the cloud. Hydrodynamic and radiative instabilities cause the cold shell to fragment as it progresses toward the center of the cloud. Similar results are seen in all of our cooling-dominated runs. As highlighted in Mellema, Kurk, & Röttgering (2002), the fragmentation is a unique result of hydrodynamic simulations of shock-cloud interactions in the cooling-dominated regime and is not observed in simulations that ignore radiative cooling (Klein, McKee, & Colella, 1994). In our simulations, the final densities of the fragments are often 3-4 orders of magnitude higher than the initial cloud density. However, as mentioned above, these final densities are generally limited by numerical resolution and are likely to underestimate the true final densities in these fragments. We are also unable to clearly resolve any interactions among either the main colliding shocks or the cold dense filaments of gas at late times.

### 4.1 Cooling Efficiency

We attempt to quantify the efficiency of the cooling processes in each of the runs presented. In order to accomplish this, we track the gas inside the cloud using a tracer fluid which is passively advected in the same manner as the density. Throughout each run the tracer distribution reflects the distribution of original cloud material. More importantly, it allows us to quantify how much of the initial cloud material cools below certain cutoff temperatures. Table 3 gives the percentages of cloud material that cool below $`T=1000`$ K and $`T=100`$ K for the equilibrium cooling curve models. The thermally unstable regime of cooling ends between $`100`$ and $`1000`$ K, so the fraction of gas with $`T<1000`$ K indicates the relative importance of cooling to the cloud evolution. Extensive star formation is unlikely, however, unless the gas is able to cool to $`T<100`$ K. The ability of the gas to reach such low temperatures in our equilibrium cooling curve models is, however, affected by numerical resolution and the neglection of self-gravity, both of which affect the peak density, and hence the cooling efficiency of the gas. Therefore, the percentage of gas that cools to below 1000 K gives a strong upper limit to the percentage that might form stars, while, when well resolved, the amount that cools to below 100 K gives a more accurate measure. These results are also represented graphically in Figure 3. We see that the equilibrium cooling curve results agree quite well with the expected transition from the non-cooling to the cooling-dominated regime, shown by the dashed line.

We also note from the results that the cooling process is generally extremely efficient throughout the cloud. In the most extreme case considered using the equilibrium cooling-curve (model E2), 77% of the gas in the cloud or $`7.7\times 10^5M_{}`$ cools below 100 K. Large fractions of cool gas were seen in all of the other models within region II. The smallest value (21%) occurs near the boundary between regions II and V, while a more typical average seen for the models is approximately 50%. The primary exception is in model E6b, in which three clouds were modeled (see below).

Table 4 reports these same percentages for the non-equilibrium chemistry models. As expected, none of these models cool below 100 K (see §3.1.2), although many of them show substantial cooling below 1000 K. Cooling is also less efficient in these models than with the equilibrium cooling-curve. In the non-equilibrium chemistry models, cooling below $`10^4`$ K is due solely to H2 emission. While H2 is an efficient coolant, the small fractions formed via gas-phase reactions (Table 4) are insufficient to lead to cooling as extensive as that due to metals, which dominate the equilibrium cooling curve models. Based upon our numerical results, it appears the cooling timescales for the non-equilibrium chemistry models are roughly a factor of 10 longer than for the equilibrium cooling curve models using the same parameters. The decreased cooling efficiency leads to a shift in the boundary between regions II and V towards higher densities and lower shock speeds, as can be seen in Figure 5. Nevertheless, large fractions of the gas contained within the clouds still cool to below 1000 K. Given higher grid resolution, self-gravity, and an extended chemical network including metals, each of which would be expected to enhance the cloud densities and therefore the H2 formation and cooling, a substantial fraction of gas might also be expected to cool below 100 K. The results indicate, therefore, that star formation may be extremely efficient in shocked clouds for which the cooling timescales are sufficiently short.

#### 4.1.1 Self-Gravity

The initial free-fall timescales for our models vary from $`7.3\times 10^6`$ Myr for $`n_{cl}=50`$ cm-3 (run E1) to $`1.6\times 10^8`$ Myr for $`n_{cl}=0.1`$ cm-3 (runs E8-E10). As explained above, our parameters were deliberately chosen such that the clouds would not be gravitationally unstable initially and that the free-fall timescale would be longer than the dynamical timescale (taken as the cloud compression time $`t_{cc}`$). Nevertheless, as we have seen, the density of the gas increases substantially in the radiative shell behind the compression shock. Therefore the local free-fall time of this compressed gas will be much shorter than its initial value. This, of course, is the basis of the shock-induced star-formation model.

For a fairly typical density increase of $`10^4`$, the local free-fall timescale will decrease by $`10^2`$. In the coldest gas, the temperature decreases by more than two orders of magnitude from its original value, and so the cold, dense regions have Jeans masses reduced by five orders of magnitude relative to the original cloud. Thus, for the higher density runs (E1-E4), the self-gravity of the gas becomes dynamically important in the dense fragments behind the compression shock. Although the Cosmos code used in this study includes an option for solving self-gravity, most of our runs did not utilize it. In the two-dimensional limit considered in the current study, we would be solving the gravity of an infinite cylindrical cloud, which is notably different than the self-gravity of a spheroidal cloud. Nevertheless, we did conduct one run (E1b) with the self-gravity option turned on. The enhanced compression due to self-gravity did allow the gas to cool more efficiently than the corresponding run without self-gravity (E1a). As a result, 69% of the gas in model E1b cooled to below 100 K, as compared to 57% in model E1a.

#### 4.1.2 Multiple Clouds

For most of this study we consider the interaction of a planar shock with a single cloud. However, our general notion of jet-induced star formation envisions the interaction of such a shock with a system of clouds within an inhomogeneous background. To highlight the possible effects of the interaction of multiple clouds, we have considered one run with a system of four clouds: two aligned perpendicular to the shock front, and two aligned parallel to it. This set up allows us to achieve higher spatial resolution by imposing symmetry boundary conditions and computing only half of the problem. Each individual cloud has the same parameters as the single cloud in model E6a ($`R_{cl}=100`$ pc, $`T_{cl}=10^4`$ K, and $`n_{cl}=1`$ cm-3). The initial cloud configuration is shown in the top of Figure 6. The size of the computational grid and the number of zones are scaled equally to maintain the same $`0.5\times 0.5`$ pc per zone resolution as our other runs.

The interaction of these clouds results in less efficient cooling throughout all of the clouds, but particularly for those in the downwind direction. The leading cloud, in fact, cools almost identically to the isolated cloud case. The discrepancy in the amount of gas found below our 100 K threshold (9% for the front cloud in run E6b versus 21% for the isolated cloud in run E6a) disappears if we change the lower threshold to 200 K. The amount of gas below that threshold is 28% for both models. This indicates that there is a significant amount of gas near 100 K in both models, but in the multi-cloud model most of it is just above this threshold.

The reduced amount of cold gas in the downwind clouds is likely due to a combination of effects. First, geometric screening of the primary shock results in less symmetric, and generally less efficient, compression for the downwind clouds. The close proximity of the clouds in this simulation also subjects them to multiple reflected shocks from their neighbors which can reheat previously cooled gas. This may be enhanced in the present work by the symmetry imposed by the two-dimensional geometry. In three dimensions, these reflected shocks would presumably be less focused. Finally, consistent with Poludnenko, Frank, & Blackman (2002), we find that the channel between the clouds acts as a de Lavalle nozel and accelerates the background gas. Since both the Rayleigh-Taylor and Kelvin-Helmholtz instabilities scale as the inverse of the background velocity, this acceleration reduces the destruction time for the downwind clouds. This effect may be reduced in three dimensions.

The results of this multi-cloud run suggest that the efficiencies we quote for our single cloud runs should generally be viewed as applicable to relatively isolated clouds. We note, however, that although reduced by a factor of a few relative to the single-cloud models, the fraction of cold gas is still significant in the model with multiple clouds. In future work we plan to present a more detailed study of the interactions of shocks with systems of radiating clouds.

### 4.2 H2 Formation

Table 4 also lists the local peak mass fraction of H2 \[max($`\rho _{H_2}/\rho `$)\] and the total H2 mass fraction ($`M_{H_2}/M_{cl}`$) for the non-equilibrium chemistry models, where $`M_{H_2}`$ is the total mass of H2 at the end of each run ($`t=1.25t_{cc}`$) and $`M_{cl}`$ is the initial cloud mass. The observed H2 mass fractions are consistent with the results of Anninos et al. (1997). Figure 7 illustrates how the H2 distribution, temperature, density, and pressure trace each other for model C1a.

#### 4.2.1 Chemistry Parameters

Along with the six combinations of $`n_{cl}`$ and $`v_{sh,b}`$ studied with the non-equilibrium chemistry model, we also explore the dependence of our results on our treatments of photoionization and photodissociation. Model C1a is our base model, comparable to the equilibrium-cooling model E3. Models C1a and C2-C6 include photoionization, photodissociation, and integrated optical depth calculations. Model C1b ignores photoionization and photodissociation. In the absence of these processes, there is also no need to calculate any optical depths. This run results in a slightly higher percentage of gas cooling below 1000 K, relative to model C1a, due to a higher mass fraction of H2. Model C1c includes photoionization and photodissociation but does not include optical depth calculations. The omission of self-shielding results in a dramatic reduction in the amount of gas which is able to cool below 1000 K and in the amount of H2 that forms, illustrating the importance of self-shielding. Model C1d increases the strength of the external photoionization field by a factor of 10. This produces a slight increase in the percentage of gas able to cool below 1000 K, due to the increased ionization fraction, which, in turn, enhances the formation of H2 via H-. Model C1e increases the H2 photoreactive destruction rate coefficient ($`k27`$) by an order of magnitude. This results in a reduction of about 14% in the amount of gas cooling below 1000 K and a dramatic reduction in the final quantity of H2.

## 5 Implications for Cloud and Jet Evolution

Our numerical studies have concentrated primarily on following the evolution of a single dense cloud being overrun by a planar shock traveling through a low density background. This picture is intended to represent what is happening on a relatively small scale ($`100`$ pc) within a much larger region of interaction ($`1`$ kpc) between a jet-induced shock and an inhomogeneous medium. We now generalize our results to predict some properties of jet-induced star-forming regions.

### 5.1 Star Formation Rates

We first wish to consider the star-formation rate for this process. To arrive at this, we need to estimate the rate at which the mass of gas in dense clumps is swept over by the shock. This is given by

$$\dot{M}_{cl}=f_{cl}\rho _{cl}v_{sh}A,$$
(24)
where $`f_{cl}`$ is the volume filling factor of the clumps and $`A`$ is the surface area of the shock front. The star-formation rate is the mass rate multiplied by the star-formation efficiency. Here it will help to consider a specific object, so we take Minkowski’s Object. From our simulations in §2.1, we estimate the shock speed inside the cloud to be $`3\times 10^8`$ cm s-1. From van Breugel et al. (1985) we take $`A20`$ kpc2, where we assume, for simplicity, that the cross-sectional area is constant. Then, for a single dense phase of $`\rho _{cl}10^{23}`$ g cm-3, we can match the observed star-formation rate of $`0.3M_{}`$ yr-1 by assuming a volume filling factor of $`f_{cl}0.03`$ and a star-formation efficiency of 0.1%.

### 5.2 Acceleration of Clouds and Mass Loading

We have already estimated the acceleration timescale of these clouds in §2.2. From this we can estimate their velocities as a function of time. In terms of the typical parameters considered in this study, we get

$$\left(\frac{v_{cl}}{\mathrm{km}/\mathrm{s}}\right)20\left(\frac{\chi }{10^3}\right)^1\left(\frac{R_{cl}}{100\mathrm{p}\mathrm{c}}\right)^1\left(\frac{v_{sh,b}}{10^3\mathrm{km}/\mathrm{s}}\right)^2\left(\frac{t}{\mathrm{Myr}}\right).$$
(25)
This value seems to overestimate the cloud velocity by a factor of a few based upon our numerical results. However, this is reasonable since the estimate assumes a constant cross-sectional area for the cloud. As we have seen, the cross-sectional area can be reduced significantly during compression. The final conclusion is that, for the parameter space studied, these clouds accelerate very slowly. We note that, for Minkowski’s Object, this is in good agreement with observations, which show very little evidence for velocity gradients in the ionized gas downstream from the jet-cloud collision region.

Along with the potential of becoming starburst regions, the nearly stationary cloud fragments may also become significant sources of mass loading for the post-shock flow (Hartquist & Dyson, 1988). Mass loading is the feeding of material into the flow and could have a significant impact on its properties, including possibly causing it to transition to the transonic regime. However, it is difficult to estimate the significance of mass loading from our simulations. Over the evolution time of our models ($`1.25t_{cc}`$), very little mass is lost from the clouds (typically $``$ 1%). This amount would certainly increase if longer evolution times were considered, but as we pointed out in §3 we can not reliably follow the evolution of these clouds beyond the current limits since our numerical resolution prevents further compression of the clouds. The question of mass loading is one we will return to in future studies.

## 6 Conclusions

We have performed two-dimensional simulations of the evolution of radiatively cooling clouds subjected to strong shocks, as might arise from galactic jets. The results of the models are summarized in Figures 3 and 5, which show the locations of the models in cloud density-shock velocity parameter space. Figure 3 presents the results for models which include equilibrium cooling with solar abundances, while Figure 5 shows the results for primordial clouds, which include non-equilibrium chemistry in which cooling at low temperatures is dominated by H2 emission. The numbers by each model indicate the percentages of the original gas that cools to below either 1000 K or 100 K. As can be seen from the figures, large fractions ($`1/3`$) of the clouds are able to cool below 1000 K. Most of this cold gas ends up in very dense filaments, which have long dynamical lifetimes. The final densities of these filaments are high enough in many cases that they are gravitationally unstable. The subsequent gravitational collapse of these regions will likely result in a burst of star formation. Our numerical results are consistent with previous work (e.g. Mellema, Kurk, & Röttgering, 2002) and have extended it by including new physics (chemistry and self-gravity) and exploring more of the parameter space of interest.

We applied our results to the specific example of Minkowski’s Object. A number of important conclusions can be drawn. First, its peculiar morphology - bright star forming region orthogonal to the jet, and fainter filamentary features downstream from there - can be easily reproduced by our simulations (Fig 1). Second, the modest amount star formation required - $`0.3M_{}`$ yr-1 for the entire object, is also easily achieved for the plausible parameter space explored by our simulations. Third, and most interestingly, we conclude that the star formation in Minkowski’s Object could be induced by a moderate velocity jet ($`9\times 10^4`$ km s-1) interacting with a collection of slightly overdense ($`10`$ cm-3), warm ($`10^4`$ K) clouds, i.e. it is NOT necessary to assume that this was an accidental collision between a jet and a preexisting gas-rich galaxy. This also suggests that the neutral hydrogen associated with Minkowski’s Object ($`3\times 10^8M_{}`$; W. van Breugel & J. van Gorkom 2003, private communication) may have cooled from the warm gas phase as a result of the radiative cooling triggered by the radio jet. An update on the observations of Minkowski’s Object will be presented in a forthcoming paper (S. D. Croft et al., in preparation).

Our results show that jets can drive radiative shocks in overdense regions in the intergalactic medium, resulting in star formation far away from the galaxies in which the jets originate. This could be an important factor in understanding the feedback of AGN on their environment. Our simulations may also be applied to radiative shocks and resulting star formation triggered by other mechanisms, such as in the wake of supernova blast waves (Preibisch et al., 2002), collisions of interstellar or intergalactic clouds (Smith, 1980), and ram pressure on gas-rich galaxies that fall into galaxy cluster atmospheres (e.g. Gavazzi et al., 2001).

Our models can be used to infer the importance of shock-induced star formation in other regimes by simply adjusting the criteria used to develop Figure 3. Our models have confirmed the reality of the division between regions II and V (cooling and non-cooling) derived using approximate analytical methods, while the other region boundaries are set by physical limits for the clouds. In future work, we shall extend our numerical simulations into regimes appropriate for other possible shock-cloud interaction scenarios.

This work was performed under the auspices of the U.S. Department of Energy by University of California, Lawrence Livermore National Laboratory under Contract W-7405-Eng-48. W.v.B. also acknowledges NASA grants GO 9779 and GO3-4150X in support of high-redshift radio galaxy research with HST and Chandra.


format: text

text:
# 1 A typical Feynman diagram for 𝐽/𝜓→𝛾ϕϕ, other diagrams are obtained by changing attachments of gluons to different places. Replacing 𝑐𝑐̄𝑔 and 𝑐𝑐̄𝛾 coupling with the corresponding dipole interaction one obtains CP violating amplitude.

CP Test in $`J/\psi \gamma \varphi \varphi `$ Decay

J.P. Ma, R.G. Ping and B.S. Zou

CCAST(World Lab.), P.O.Box 8730, Beijing 100080, China;
Institute of High Energy Physics, Academia Sinica, Beijing 100039, China;
Institute of Theoretical Physics , Academia Sinica, Beijing 100080, China

## Abstract

We propose to test CP symmetry in the decay $`J/\psi \gamma \varphi \varphi `$, for which large data sample exists at BESII, and a data sample of $`10^{10}`$ $`J/\psi `$’s will be collected with BESIII and CLEO-C program. We suggest some CP asymmetries in this decay mode for CP test. Assuming that CP violation is introduced by the electric- and chromo-dipole moment of charm quark, these CP asymmetries can be predicted by using valence quark models. Our work shows a possible way to get information about the electric- and chromo-dipole moment of charm quark, which is little known. Our results show that with the current data sample of $`J/\psi `$, electric- and chromo-dipole moment can be probed at order of $`10^{13}ecm`$. In the near future with a $`10^{10}`$ data sample, these moments can be probed at order of $`10^{14}ecm`$.

PACS numbers: 11.30.Er,13.25.Gv

CP violation was first discovered in $`K`$-system in the last century. At the beginning of this century, CP violation was also found in $`B`$-system thanks to two $`B`$-factories. Except these, there is no experimental indication for CP violation. With our current knowledge interactions between elementary particles are described by the Standard Model(SM). In SM CP symmetry is violated by complex CKM matrix elements. Although SM is successful and can explain CP violation effects found in $`K`$\- and $`B`$-system, it is generally believed that SM is only an effective theory for energy scale around $`100`$GeV. It is expected that new physics beyond SM will emerge. Therefore, it is important to find CP violation somewhere else, where SM predicts no or tiny CP violation, because this can be an indication for new physics.

It is expected that effect of CP violation is very small in general. Hence, Experimental study of CP violation needs large data samples. With progress at BEPC and CLEO-C large data samples of $`J/\psi `$ decays exist or will be collected. At BEPC with BESII detector there are already $`5.8\times 10^7`$ $`J/\psi `$-decay events. With upgraded BEPC and CLEO-C program, a data sample of $`10^910^{10}`$ $`J/\psi `$ will be collected. These huge data samples are very suited for CP test. But, not every decay mode of $`J/\psi `$ can be used for CP test. For $`J/\psi `$ decay into a particle and its antiparticle, CP test is not possible if these particles are spinless or their polarizations are not observed. CP test is only possible if polarizations of decay products are measured. The decay $`J/\psi \mathrm{\Lambda }\overline{\mathrm{\Lambda }}`$ is an example, where the polarizations can be measured though subsequential decays of $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$. CP test with this decay mode has been studied in . CP test can be made for three-body decays, even without knowing polarizations of decay products. In this letter we propose to study CP test in the decay $`J/\psi \gamma \varphi \varphi `$. This decay can provide useful information about electric dipole moment (EDM) of charm quark. The reason for choosing this channel is because $`\varphi `$ is a very narrow resonance just above $`K\overline{K}`$ threshold and can be clearly identified by its $`K^+K^{}`$ decay mode in experiment. In principle $`J/\psi \gamma \rho \rho `$ could also serve the purpose, but experimentally the broad width of $`\rho `$ meson makes it impossible to get clean sample for the channel. By using large charmonia data samples it has been also discussed in  about possible CP test, e.g., $`\psi ^{}J/\psi \pi \pi `$, and in  about parity test in the decay mode $`J/\psi \gamma \mathrm{}^+\mathrm{}^{}`$.

From SM EDM of quarks and leptons are very small(See Reviews in , and references therein). If EDM of quark is found to be nonzero, it is likely an indication for new physics. Because the operator for EDM does not converse helicities of quarks, its effect is suppressed in a high energy process by a factor $`m_q/E`$, where $`m_q`$ is the quark mass and $`E`$ is a large energy scale. For light quarks, useful information can be obtained through measurement of EDM of neutron . So far there is no experimental information about EDM’s of heavy quarks, like charm- and bottom quark. $`J/\psi `$ decays can provide information of EDM of charm quark and has the advantage that the effect of EDM will be not suppressed, because the large energy scale is around $`m_c`$. Since in the decay a $`c\overline{c}`$ pair is annihilated into a photon and gluons, it also provides a way to detecting chromodipole moment of charm quark. In this letter we assume that CP violation is introduced by electric- and chromo-dipole moment and study their effect in the decay. The moments are defined by effective Lagrangian:

$`L_{CP}=i{\displaystyle \frac{d_c}{2}}\overline{c}\gamma _5\sigma _{\mu \nu }F^{\mu \nu }ci{\displaystyle \frac{\stackrel{~}{d}_c}{2}}\overline{c}\gamma _5\sigma _{\mu \nu }G^{\mu \nu }c,`$ (1)
where $`d_c`$ is the electric dipole moment, $`\stackrel{~}{d}_c`$ is the chromodipole moment. With the decay mode, it is also possible to detect the electric dipole moment of $`s`$-quark. But its effect is suppressed by $`m_s/m_c`$. In this letter we will only consider dipole moments of $`c`$-quark.

Unfortunately, predictions for the decay $`J/\psi \gamma \varphi \varphi `$ can not be made without models, especially, because there is no reliable method in theory how to determine the formation of bound states like $`\varphi `$. For $`J/\psi `$ a good approximation is to use a nonrelativistic wave-function for describing $`J/\psi `$ as a bound state of a $`c\overline{c}`$ pair, because $`c`$-quark can be taken as a heavy quark. In this work we will also take this approximation for $`\varphi `$ meson by taking it as a nonrelativistic bound state of a $`s\overline{s}`$ quark pair, where the $`s`$-quark should be taken as a constituent quark, but not as a current quark. The constituent $`s`$-quark mass should be taken as $`m_sm_\varphi /2`$. This approximation was employed in radiative decay into a hadron. In this approximation, the decay amplitude is factorized into a partonic part which is for the process $`c\overline{c}\gamma s\overline{s}s\overline{s}`$, and a nonperturbative part which is related to wave functions of hadrons at the origin. The partonic part can be calculated with perturbative theory. The wave functions at the origin can be determined by experiment. In this letter we will take this approximation. It should be noted that with this approximation reliable predictions for various distributions may not be obtained. But one may expect that for our integrated observables defined in below it can be a good approximation. Especially, our integrated observables will not depend on wave functions at the origin.

We consider the decay in the rest-frame of $`J/\psi `$

$$J/\psi (P)\gamma (k)+\varphi (p_1)+\varphi (p_2),$$
(2)
where momenta are given in brackets. Because the two $`\varphi `$ mesons are identical particles, we require $`p_1^0>p_2^0`$ to distinguish them in experiment. The decay amplitude can be written as:

$$𝒯(J/\psi \gamma \varphi \varphi )=A_i(𝐩_\mathrm{𝟏},𝐩_\mathrm{𝟐})\epsilon _i,$$
(3)
where $`\epsilon _i(i=1,2,3)`$ is the spacial component of the polarization vector of $`J/\psi `$. With the amplitude, a density matrix for the decay can be defined as:

$$R_{ij}(𝐩_\mathrm{𝟏},𝐩_\mathrm{𝟐})=\underset{spin}{}A_i(𝐩_\mathrm{𝟏},𝐩_\mathrm{𝟐})A_j^{}(𝐩_\mathrm{𝟏},𝐩_\mathrm{𝟐}).$$
(4)
The sum is over spins of all particles in the final state. This matrix can be decomposed in terms of tensor built with $`𝐩_\mathrm{𝟏}`$, $`𝐩_\mathrm{𝟐}`$, $`\delta _{ij}`$ and $`\epsilon _{ijk}`$. A similar decomposition can be found in . The density matrix contains all information about the decay. Any experimental observable can be predicted if we know the density matrix for production. We consider the situation where $`J/\psi `$ is produced through $`e^+e^{}`$ annihilation in their center-mass frame, as at BEPC or CLEO-C. The density matrix $`\rho _{ij}`$ for the production takes the form:

$$\rho _{ij}=\frac{1}{3}\delta _{ij}\frac{1}{2}(\widehat{k}_{+i}\widehat{k}_{+j}\frac{1}{3}\delta _{ij}),$$
(5)
where $`\widehat{k}_+`$ denotes the moving direction of $`e^+`$. With these density matrices, the expectation value of any experimental observable $`O`$ constructed with $`𝐤_+`$, $`𝐩_\mathrm{𝟏}`$ and $`𝐩_\mathrm{𝟐}`$ can be predicted by

$$O=\frac{1}{𝒩}𝑑\mathrm{\Gamma }(2\pi )^4\delta (Pp_1p_2k)O\rho _{ij}R_{ji},$$
(6)
where $`d\mathrm{\Gamma }`$ is the differential phase-space of $`\gamma \varphi \varphi `$ with $`p_1^0>p_2^0`$, $`𝒩`$ is a normalization factor so that $`1=1`$.

Symmetries will constrain the form of the density matrix. If CP symmetry holds, we have

$$R_{ij}(𝐩_\mathrm{𝟏},𝐩_\mathrm{𝟐})=R_{ij}(𝐩_\mathrm{𝟏},𝐩_\mathrm{𝟐}).$$
(7)
For CP test one can check this relation to see if it holds. A convenient way to check this is to use some integrated observables. We find two CP odd observables which we are particularly interested in:

$$O_1=\widehat{𝐤}_+\widehat{𝐩}_\mathrm{𝟏}\widehat{𝐤}_+(\widehat{𝐩}_\mathrm{𝟏}\times \widehat{𝐩}_\mathrm{𝟐}),O_2=\widehat{𝐤}_+\widehat{𝐩}_\mathrm{𝟐}\widehat{𝐤}_+(\widehat{𝐩}_\mathrm{𝟏}\times \widehat{𝐩}_\mathrm{𝟐}),$$
(8)
where $`\widehat{𝐩}_\mathrm{𝟏}`$ or $`\widehat{𝐩}_\mathrm{𝟐}`$ are the direction of $`𝐩_\mathrm{𝟏}`$ or $`𝐩_\mathrm{𝟐}`$ respectively. From these oberservables, One can define CP -asymmetry as

$$B_i=\theta (O_i)\theta (O_i)(i=1,2),$$
(9)
where $`\theta (x)=1`$ if $`x>0`$ and is zero if $`x<0`$. If these asymmetries are not zero, CP symmetry is violated. In this letter, we will make predictions for these asymmetries by taking CP violating interactions as given in Eq.(1).

As discussed before, we will use nonrelativistic wave-functions for bound states. A typical Feynman diagram is given in Fig.1. From which one can read out the contribution to the scattering amplitude for $`c(q_1)+\overline{c}(q_2)\gamma (k)+s(k_1)+\overline{s}(k_2)+s(k_3)+\overline{s}(k_4)`$, which takes the form:

$$𝒯_1(c(q_1)+\overline{c}(q_2)\gamma (k)+s(k_1)+\overline{s}(k_2)+s(k_3)+\overline{s}(k_4))=\overline{v}(q_2)C^{ab,\mu \nu }u(q_1)\overline{u}(k_1)\gamma _\mu T^av(k_4)\overline{u}(k_3)\gamma _\nu T^bv(k_2),$$
(10)
where $`C^{ab,\mu \nu }`$ is a matrix with Dirac- and color indices. To obtain the contribution from Fig.1 to the hadronic decay amplitude, one needs to project the quark pairs into hadron states. The projection is made by the replacement:

$`u_i(q_1)\overline{v}_j(q_2)`$ $``$ $`{\displaystyle \frac{1}{3(2m_c)^{3/2}}}\psi _{J/\psi }(0)\left[(\gamma q_1+m_c)\gamma \epsilon (P)(\gamma q_2+m_c)\right]_{ij},`$
$`v_i(k_2)\overline{u}_j(k_1)`$ $``$ $`{\displaystyle \frac{1}{3(2m_s)^{3/2}}}\psi _\varphi ^{}(0)\left[(\gamma k_2+m_s)\gamma \epsilon ^{}(p_1)(\gamma k_1+m_s)\right]_{ij},`$
$`v_i(k_4)\overline{u}_j(k_3)`$ $``$ $`{\displaystyle \frac{1}{3(2m_s)^{3/2}}}\psi _\varphi ^{}(0)\left[(\gamma k_4+m_s)\gamma \epsilon ^{}(p_2)(\gamma k_3+m_s)\right]_{ij},`$ (11)
where $`ij`$ stand for Dirac- and color indices, $`\epsilon (P)`$, $`\epsilon (p_1)`$ and $`\epsilon (p_2)`$ is the polarization vector of $`J/\psi (P)`$, $`\varphi (p_1)`$ and $`\varphi (p_2)`$, respectively. $`\psi _{J/\psi }`$($`\psi _\varphi `$) is the wave function of $`J/\mathrm{\Psi }`$($`\varphi `$). In the projection momenta of quarks and masses should be taken as:

$`q_1`$ $`=`$ $`q_2={\displaystyle \frac{P}{2}},m_c={\displaystyle \frac{M_{J/\psi }}{2}},`$
$`k_1`$ $`=`$ $`k_2={\displaystyle \frac{p_1}{2}},m_s={\displaystyle \frac{M_\varphi }{2}},`$
$`k_3`$ $`=`$ $`k_4={\displaystyle \frac{p_2}{2}}.`$ (12)
It should be noted that from the amplitude in Eq.(10) there is another contribution to the hadronic amplitude which is obtained by exchanging $`p_1p_2`$. The calculation is straightforward. The total decay amplitude is too complicated to be given here. A numerical program for the amplitude can be obtained from authors.

With the approximation explained in detail, numerical results can be obtained. We use $`\alpha _s(m_c)=0.28`$ and masses which are fixed by experimental values of hadrons given in Eq.(12). We obtain the decay width:

$$\mathrm{\Gamma }(J/\psi \gamma \varphi \varphi )=1.83\times 10^3|\psi _{J/\psi }(0)|^2|\psi _\varphi (0)|^4\mathrm{GeV}^8.$$
(13)
If we use the value of $`|\psi _{J/\psi }(0)|^2=0.055\mathrm{GeV}^3`$, determined from study with potential models and lattice QCD, and experimental results of $`\mathrm{\Gamma }_{tot}=87\mathrm{K}\mathrm{e}\mathrm{V}`$ and $`Br(J/\psi \gamma \varphi \varphi )=4.0\times 10^4`$, we can obtain

$$|\psi _\varphi (0)|^20.019\mathrm{GeV}^3.$$
(14)
It should be noted that in our approximation the CP asymmetries in Eq.(9) do not depend on these wave functions at the origin. With the same parameters we obtain:

$`B_1`$ $`=`$ $`4.2\left[{\displaystyle \frac{d_c}{10^{10}ecm}}\right]1.2\left[{\displaystyle \frac{\stackrel{~}{d}_c}{10^{10}ecm}}\right],`$
$`B_2`$ $`=`$ $`3.9\left[{\displaystyle \frac{d_c}{10^{10}ecm}}\right]+1.3\left[{\displaystyle \frac{\stackrel{~}{d}_c}{10^{10}ecm}}\right].`$ (15)
A statistic sensitivity to $`d_c`$ and $`\stackrel{~}{d}_c`$ can be determined from these results by requiring that the asymmetry generated by these dipole moments should be larger than the statistical error. The statistical error for these asymmetries is given by

$$\delta B_i=\sqrt{\frac{1}{N_{events}}},$$
(16)
where $`N_{events}`$ is the number of the decay events. With the $`5.8\times 10^7J/\psi `$ data sample at BES, the sensitivities of our CP asymmetries to these dipole moments are

$$d_c1.4\times 10^{13}ecm,\stackrel{~}{d}_c4.5\times 10^{13}ecm.$$
(17)
With a $`10^{10}`$ data sample which will be collected in the near future, the sensitivities are:

$$d_c1.2\times 10^{14}ecm,\stackrel{~}{d}_c3.6\times 10^{14}ecm.$$
(18)
It is instructive to compare these sensitivities to the upper bound of EDM of neutron and electron from experiment. The upper bounds read:

$`d_n`$ $`<`$ $`6.3\times 10^{26}ecm,`$
$`d_e`$ $`<`$ $`4.3\times 10^{27}ecm.`$ (19)
The implication of these upper bounds for physics beyond SM has been studied extensively(See e.g.,  and references therein). Predictions from any reasonable extension of SM for EDM of light quarks and electron should not be larger than these upper bounds. Comparing these upper bounds, the sensitivities are too large. However, in some extensions of SM predicted EDM of a quark is proportional to the mass of the quark(e.g., see ). It is possible to obtain EDM of $`c`$-quark which is comparable to these sensitivities. Hence, CP test proposed in this letter will provide information for possible new physics beyond SM.

To summarize: In this letter we have proposed to test CP symmetry in the decay $`J/\psi \gamma \varphi \varphi `$ by using large data samples of $`J/\psi `$ which is already collected at BES and will be collected with BESIII and CLEO-C program. We introduce some CP asymmetries for CP test with this decay mode. Assuming that CP violation is introduced by the electric- and chromo-dipole moment of charm quark, we give predictions for these asymmetries by using wave functions for hadrons. Experimental information is available for electric dipole moments of light quarks through measurement of the electric dipole moment of neutron, but there is no information for electric dipole moments of heavy quarks. Our work shows a possible way to obtain this information. With the current data sample of $`J/\psi `$ at BES, the electric- and chromo-dipole moment can be probed at order of $`10^{13}ecm`$. In the near future with a $`10^{10}`$ data sample, these moments can be probed at order of $`10^{14}ecm`$.

Acknowledgements

One of us(J.P. Ma) would like to thank Prof. X.G. He and O. Nachtmann for useful discussions. This work is supported by National Nature Science Foundation of P. R. China.


text:
# Superconducting qubit storage and entanglement with nanomechanical resonators

## Abstract

We describe a quantum computational architecture based on integrating nanomechanical resonators with Josephson junction phase qubits, with which we implement single- and multi-qubit operations. The nanomechanical resonator is a GHz-frequency, high-quality-factor dilatational resonator, coupled to the Josephson phase through a piezoelectric interaction. This system is analogous to one or more few-level atoms (the Josephson qubits) in a tunable electromagnetic cavity (the nanomechanical resonator). Our architecture combines the best features of solid-state and cavity-QED approaches, and may make possible multi-qubit processing in a scalable, solid-state environment.

The lack of scalable qubit architectures, with sufficiently long quantum-coherence lifetimes and a suitably controllable entanglement scheme, remains the principal roadblock to building a large-scale quantum computer. Superconducting devices exhibit robust macroscopic quantum behavior Makhlin et al. (2001). Recently, there have been exciting demonstrations of long-lived Rabi oscillations in current-biased Josephson junctions Yu et al. (2002); Martinis et al. (2002), subsequently combined with a two-qubit coupling scheme Berkley et al. (1998), and in parallel, demonstrations of Rabi oscillations and Ramsey fringes in a Cooper-pair box Nakamura et al. (1999, 2002); Vion et al. (2002). These accomplishments have generated significant interest in the potential for Josephson-junction-based quantum computation Leggett (2002). Coherence times $`\tau _\phi `$ up to 5$`\mu `$s have been reported in the current-biased devices Yu et al. (2002), with corresponding quality factors $`Q_\phi \tau _\phi \mathrm{\Delta }E/h`$ of the order of $`10^5`$, yielding sufficient coherence to perform many logical operations. Here $`\mathrm{\Delta }E`$ is the qubit energy level spacing.

In this paper, we describe an architecture in which ultrahigh-frequency resonators coherently couple two or more current-biased Josephson junctions, where the superconducting “phase qubits” are formed from the energy eigenstates of the junctions. We show that the system is analogous to one or more few-level atoms (the Josephson junctions) in a tunable electromagnetic cavity (the resonator), except that here we can individually tune the energy level spacing of each atom, and control the electromagnetic interaction strength.

Other investigators have proposed the use of electromagnetic Shnirman et al. (1997); Makhlin et al. (1999); Mooij et al. (1999); Makhlin et al. (2000); You et al. (2002); Plastina and Falci (2003); Blais et al. (2003); Smirnov and Zagoskin  or superconducting Buisson and Hekking (2001); Marquardt and Bruder (2001) resonators to couple Josephson junctions together. The use of nanomechanical resonators to mediate multi-qubit operations has not to our knowledge been described previously, although an approach to create entangled states of a single resonator has been proposed Armour et al. (2002). The use of mechanical as opposed to electromagnetic resonators has the advantage that potentially much higher quality factors can be achieved Yang et al. (2000), with significantly smaller dimensions, enabling a truly scalable approach.

Our implementation uses large-area current-biased Josephson junctions, with capacitance $`C`$ and critical current $`I_0`$; a circuit model is shown in Fig. 1. The largest relevant energy is the Josephson energy $`E_\mathrm{J}\mathrm{}I_0/2e`$, with a charging energy $`E_\mathrm{c}(2e)^2/2CE_\mathrm{J}`$. The dynamics of the Josephson phase difference $`\delta `$ is that of a particle of mass $`M=\mathrm{}^2C/4e^2`$ moving in an effective potential $`U(\delta )E_\mathrm{J}(\mathrm{cos}\delta +s\delta )`$, for bias current $`s=I_\mathrm{b}/I_0`$ Barone and Paterno (1982); Fulton and Dunkleberger (1974). For bias currents $`s<1`$, the potential $`U(\delta )`$ has metastable minima, separated from the continuum by a barrier $`\mathrm{\Delta }UU(\delta _{\mathrm{max}})U(\delta _{\mathrm{min}})(4\sqrt{2}/3)E_\mathrm{J}(1s)^{3/2}`$ for $`s1^{}`$, as shown in Fig. 1. The curvature $`U^{\prime \prime }(\delta )`$ defines the small-amplitude plasma frequency $`\omega _\mathrm{p}\sqrt{U^{\prime \prime }(\delta _{\mathrm{min}})/M}=\omega _{\mathrm{p0}}(1s^2)^{1/4}`$, with $`\omega _{\mathrm{p0}}=\sqrt{2eI_0/\mathrm{}C}=\sqrt{2E_\mathrm{c}E_\mathrm{J}}/\mathrm{}`$. The Hamiltonian for the junction phase difference is $`H_\mathrm{J}=P^2/2M+U(\delta )`$, with $`P=i\mathrm{}d/d\delta `$ the momentum operator. The junction’s zero-voltage state corresponds to the phase “particle” trapped in one of the metastable minima.

The lowest two quasi-bound states in a local minimum, $`|0`$ and $`|1`$, define the phase qubit. State preparation is typically carried out with $`s`$ just below unity, in the range $`s=0.950.99`$, where $`U(\delta )`$ is strongly anharmonic, and for which there are only a few quasibound states Martinis et al. (2002); Berkley et al. (1998). The anharmonicity allows state preparation from a classical radiofrequency (rf) field, as then the frequency of the classical field can be set to couple to only the lowest two states. In our scheme, by contrast, single quanta are exchanged between the junction and the resonator, so anharmonicity is not necessary; we find it convenient to work with $`s`$ between 0.5 and 0.9.

We focus here on coupling a single resonator to two Josephson qubits; extensions to larger systems will be considered in later work. The two-junction circuit is shown in Fig. 2. The disk-shaped element is the nanomechanical resonator, consisting of a single-crystal piezoelectric disc sandwiched between two metal plates, and the junctions are the crossed boxes on either side of the resonator, interrogated by high-impedance circuits Martinis et al. (2002).

The phase qubit state $`|0`$ of a single junction is prepared by waiting for any excited component to decay. The pure state $`|1`$, or a superposition state $`\alpha |0+\beta |1`$, is prepared by adding a classical rf current $`I_{\mathrm{rf}}`$ to the bias, $`I_{\varphi 1}(t)=I_{\mathrm{dc}}+I_{\mathrm{rf}}^\mathrm{c}\mathrm{cos}(\omega _{\mathrm{rf}}t)+I_{\mathrm{rf}}^\mathrm{s}\mathrm{sin}(\omega _{\mathrm{rf}}t)`$. Both $`I_{\mathrm{dc}}`$ and $`I_{\mathrm{rf}}^{\mathrm{s},\mathrm{c}}`$ vary slowly compared to $`\mathrm{}/\mathrm{\Delta }E`$. When $`\omega _{\mathrm{rf}}`$ is near resonance with the level spacing $`\mathrm{\Delta }E/\mathrm{}`$, the qubit will undergo Rabi oscillations, allowing the controlled preparation of linear combinations of $`|0`$ and $`|1`$.

The nanomechanical resonator is designed with a fundamental thickness resonance frequency $`\omega _0/2\pi 110\mathrm{GHz}`$, with quality factor $`Q10^510^6`$. Piezoelectric dilatational resonators with resonance frequencies in this range, and quality factors of $`10^3`$ at room temperature, have been fabricated from sputtered AlN Ruby and Merchant (1994); Ruby et al. (2001). Single-crystal AlN can also be grown by chemical vapor deposition Cleland et al. (2001). Our simulations are based on such a resonator, with a diameter $`d=1.16\mu `$m and thickness $`b=0.5\mu m`$ pie . Such resonators can be used to coherently store a qubit state prepared in a current-biased Josephson junction, return it to that junction, or transfer it to another junction, as well as entangle two or more junctions. These operations are performed by tuning the energy level spacing $`\mathrm{\Delta }E`$ into resonance with $`\mathrm{}\omega _0`$, generating electromechanical Rabi oscillations.

Referring to Fig. 2, the total bias current of junction 1 is $`I_{\mathrm{dc1}}+I_{\mathrm{res}}`$, where $`I_{\mathrm{res}}`$ is the current through the resonator from that junction. A simple model for the resonator allows us to write $`I_{\mathrm{res}}=C_{\mathrm{res}}(\dot{V}+h_{33}b\dot{U})`$, where $`C_{\mathrm{res}}`$ is the resonator geometric capacitance, $`h_{33}`$ the relevant piezoelectric coupling constant Auld (1990), $`\dot{V}`$ the rate of voltage change, and $`\dot{U}`$ the rate of change of the mechanical strain. The current $`I_{\mathrm{res}}`$ is partly due to the capacitance $`C_{\mathrm{res}}`$ and partly due to the piezoelectrically-coupled strain $`U`$. $`C_{\mathrm{res}}`$, in parallel with the junction capacitance $`C`$, renormalizes the mass $`M`$ to $`\stackrel{~}{M}=\mathrm{}^2\stackrel{~}{C}/4e^2`$, where $`\stackrel{~}{C}=C+C_{\mathrm{res}}`$.

With the resonator coupled to the superconducting phase through the voltage $`V`$, the Hamiltonian for the combined junction-resonator system is $`H=H_\mathrm{J}+H_{\mathrm{res}}+\delta H`$. Here $`H_{\mathrm{res}}=\mathrm{}\omega _0a^{}a`$ is the Hamiltonian of the isolated resonator, where we have quantized the resonator displacement field with creation (destruction) operators $`a^{}`$ ($`a`$), and only included the fundamental dilatational mode. $`\delta H`$ is the phase-resonator interaction,

$$\delta H=\frac{\mathrm{}C_{\mathrm{res}}bh_{33}}{2e(1\eta )}\delta \dot{U}\delta =ig(aa^{})\delta ,$$
(1)
where $`\eta =0.054`$ and the coupling constant $`g`$ is

$$g=\frac{\mathrm{}^{3/2}C_{\mathrm{res}}h_{33}\sqrt{\omega _0}}{(1\eta )e\sqrt{\rho \pi bd^2/4}}.$$
(2)
For our model resonator $`g0.820\mu `$eV.

In the junction eigenstate basis, the junction Hamiltonian is $`H_\mathrm{J}=_mϵ_mc_m^{}c_m`$, with creation (destruction) operators $`c_m^{}`$ ($`c_m`$) acting on the phase qubit states. The interaction Hamiltonian is

$$\delta H=ig\underset{mm^{}}{}m|\delta |m^{}c_m^{}c_m^{}(aa^{}).$$
(3)
The eigenstates of the noninteracting Hamiltonian $`H_0=H_\mathrm{J}+H_{\mathrm{res}}`$ are $`|mn|m_\mathrm{J}|n_{\mathrm{res}}`$, with energies $`E_{mn}=ϵ_m+\mathrm{}\omega _0n`$, where $`n`$ is the resonator occupation number. An arbitrary state can be expanded as $`|\mathrm{\Psi }(t)=_{mn}c_{mn}(t)|mn\mathrm{exp}(iE_{mn}t/\mathrm{})`$.

The full Hamiltonian is equivalent to a few-level atom in an electromagnetic cavity. The cavity “photons” are phonons, which interact with the “atoms” (here the Josephson junctions) via the piezoelectric effect. This analogy allows us to adapt quantum-information protocols developed for cavity-QED to our architecture.

We first show that we can coherently transfer a qubit state from a junction to a resonator, using the adiabatic approximation combined with the rotating-wave approximation (RWA) of quantum optics Scully and Zubairy (1997). We assume that the bias current $`s`$ changes slowly on the time scale $`\mathrm{}/\mathrm{\Delta }E`$, and work at temperature $`T=0`$. The RWA is valid when $`\mathrm{\Delta }E`$ and $`\mathrm{}\omega _0`$ are close on the scale of $`\mathrm{}\omega _0/Q_{\mathrm{res}}`$, and when the interaction strength $`g\mathrm{\Delta }E`$. At time $`t=0`$, we prepare the resonator in the state $`|0_{\mathrm{res}}`$. In the RWA, neglecting relaxation, we obtain the amplitude evolution

$`i\mathrm{}_tc_{0n}`$ $`=`$ $`ig\sqrt{n}0|\delta |1e^{i\omega _\mathrm{d}t}c_{1,n1}`$
$`i\mathrm{}_tc_{1n}`$ $`=`$ $`ig\sqrt{n+1}1|\delta |0e^{i\omega _\mathrm{d}t}c_{0,n+1},`$ (4)
where $`\omega _\mathrm{d}\omega _0\mathrm{\Delta }E/\mathrm{}`$ is the resonator–qubit detuning. We integrate to find the reduced density matrices $`\rho _\mathrm{J}(t)`$ (in the qubit subspace) and $`\rho _{\mathrm{res}}(t)`$ (in the zero- and one-phonon resonator subspace). The junction phase is initially prepared in the pure state $`\alpha |0_\mathrm{J}+\beta |1_\mathrm{J}`$, corresponding to the reduced density matrix

$$\rho _\mathrm{J}(0)=\left[\begin{array}{cc}|\alpha |^2& \alpha \beta ^{}\\ \alpha ^{}\beta & |\beta |^2\end{array}\right].$$
(5)
We allow the junction and resonator to interact on resonance for a time $`\mathrm{\Delta }t=\pi /\mathrm{\Omega }_d`$, where the Rabi frequency is $`\mathrm{\Omega }_d=(\mathrm{\Omega }_0^2+\omega _\mathrm{d}^2)^{1/2}`$, in terms of the tuned (resonant) value $`\mathrm{\Omega }_0=2g|0|\delta |1|/\mathrm{}`$. After the interval $`\mathrm{\Delta }t`$, the *resonator* is found to be in the same pure state,

$$\rho _{\mathrm{res}}(\pi /\mathrm{\Omega }_d)=\left[\begin{array}{cc}|\alpha |^2& \alpha \beta ^{}e^{i\pi \omega _0/\mathrm{\Omega }_d}\\ \alpha ^{}\beta e^{i\pi \omega _0/\mathrm{\Omega }_d}& |\beta |^2\end{array}\right],$$
(6)
apart from expected phase factors. The phase qubit state has been swapped with that of the resonator. The cavity-QED analog of this operation has been demonstrated experimentally in Ref. Maître et al. (1997).

To assess the limitations of the RWA, we also numerically integrated the exact amplitude equations

$$i\mathrm{}\dot{c}_{mn}=\underset{m^{}n^{}}{}mn|\delta H|m^{}n^{}e^{i(E_{mn}E_{m^{}n^{}})t/\mathrm{}}c_{m^{}n^{}}.$$
(7)
The Josephson junction had parameters corresponding Ref. Martinis et al. (2002), $`E_\mathrm{J}=43.05\mathrm{meV}`$ and $`E_\mathrm{c}=53.33\mathrm{neV}`$. We used a 4th-order Runge-Kutta method with a time step of $`10\mathrm{fs}`$. Our main result is shown in Fig. 3. The qubit transfer depends sensitively on the shape of the profile $`s(t)`$, which starts at $`s=0.50`$, and is then adiabatically changed to the resonant value $`s=0.928`$. We find that the time during which $`s`$ changes should be at least exponentially localized. This can be understood by recalling that the RWA requires the qubit to be exactly in resonance with the resonator (in the $`Q\mathrm{}`$ limit). Therefore one must bring the system into resonance as quickly as possible without violating adiabaticity. The power-law tails associated with an arctangent function, for example, lead to large deviations from the desired behavior, shown in Fig. 3(b). The result in Fig. 3(a) was obtained using trapezoidal profiles with a cross-over time of $`0.5\mathrm{ns}`$. All quasibound junction states were included in the calculation, and convergence with the resonator’s Hilbert space dimension was obtained. The junction is held in resonance for half a Rabi period $`\pi /\mathrm{\Omega }_d`$, during which energy is exchanged at the Rabi frequency. The systems are then brought out of resonance. The final state amplitudes are given in Table 1, and are quite close to the RWA results.

To pass a qubit state $`\alpha |0+\beta |1`$ from junction 1 to junction 2, the state is loaded into the first junction and the bias current changed to bring the junction into resonance with the resonator for half a Rabi period. This writes the state $`\alpha |0+\beta |1`$ into the resonator. After the first junction is taken out of resonance, the second junction is brought into resonance for half a Rabi period, passing the state to the second junction. We have simulated this operation numerically, assuming two identical junctions coupled to the resonator described above. The results are shown in Fig. 4 and Table 2, where $`c_{m_1m_2n}`$ is the probability amplitude (in the interaction representation) to find the system in the state $`|m_1m_2n`$, with $`m_1`$ and $`m_2`$ labelling the states of the two junctions.

We can prepare an entangled state of two junctions by bringing the first junction into resonance with the resonator for 1/4th of a Rabi period Hagley et al. (1997), which, according to our RWA analysis, produces the state $`(|100|001)/\sqrt{2}`$. After bringing the second junction into resonance for half a Rabi period, the state of the resonator and second junction are swapped, leaving the system in the state $`(|100|010)/\sqrt{2}`$ with a probability of 0.987, where the resonator is in the ground state and the junctions entangled, as demonstrated in Fig. 5. Using the cavity-QED analogy, it will be possible to transfer the methodology developed for the standard two-qubit operations, in particular controlled-NOT logic, to this system, using mostly existing technology and demonstrated techniques.

Acknowledgements. ANC and MRG were supported by the Research Corporation, and MRG was by NSF CAREER Grant No. DMR-0093217.