text:
# Tidal Torques and Galactic Warps

## 1 Introduction

Many edge-on disk galaxies show integral-sign or S-shaped warps, where the majority of the disk is planar but where the outer region of the disk lies above that plane on one side of the galaxy and below the plane on the other (Binney, 1992). The Milky Way is warped both in neutral hydrogen (e.g. Diplas & Savage, 1991) and in the stellar distribution (Reed, 1996; López-Corredoira et al, 2002b). Reshetnikov & Combes (1998) estimate that half of all disk galaxies have optical warps, and most HI disks which extend beyond the optical disks appear to be warped (Bosma, 1981; Briggs, 1990; Christodoulou et al, 1993).

A number of methods have been proposed for creating and maintaining warps. Several authors (Lynden-Bell, 1965; Sparke & Casertano, 1988) have suggested that the system of discrete particles which make up the disk may have normal bending modes that could be excited. Battaner et al (1990) have investigated magnetic fields as a cause of warps, while Toomre (1983) and Dekel & Shlosman (1983) suggested that disks askew in flattened dark matter halos could develop long-lived warps, assuming the radial profile of the halo were appropriately fine-tuned.

Following on the idea that infalling material will shift the angular momentum of a galaxy (Ryden, 1988; Quinn & Binney, 1992), Ostriker & Binney (1989) studied how a disk of massive rings reacts to a slewing disk potential. They found that the reaction of the disk depends largely on its surface density, with warps appearing in regions of low surface density. Also motivated by the cosmic infall of material with differing angular momentum, Debattista & Sellwood (1999) found that when the angular momentum of a halo and disk are misaligned, dynamical friction transfers angular momentum between them in such a way as to produce a long-lived warp.

More recently, López-Corredoira et al (2002a) examined the torques produced by the transfer of angular momentum from gas falling onto a galactic disk, and the warping of the disk in response to these torques as well as the internal torques due to the interaction of different rings within the disk. They find a disk response which mirrors the Milky Way warp, though the effects of embedding the disk in a dark halo are as yet insufficiently modelled.

In a cosmological setting, a galactic disk experiences gravitational torques from three different sources: external galaxies, non-spherically-symmetric substructure in the halo, and misalignment between the disk and halo; if the dark matter halo is not spherical, as suggested both by simulations (e.g. Frenk et al, 1985; Katz, 1991; Dubinski & Carlberg, 1991; Cole & Lacey, 1996) and observations (Sackett, 1999), then it will exert a gravitational torque on a disk not in its symmetry plane.

In this paper, we evaluate the effect of the gravitational tidal torques a typical galactic disk experiences from its misalignment with the halo, and study whether these torques provide a possible origin for warped disks. We first develop a framework for how orbits in a disk react to torques. We then investigate the form and magnitude of torques due to flattened misaligned halos. The effect of these torques on massless disks is examined by analytic and numerical techniques. Finally, we perform N-body simulations of massive disks resembling the Milky Way and investigate how the self-gravity of the disk may affect its reaction to cosmological torques.

## 2 Origin of Torques on Galactic Disks

### 2.1 Framework for Torqued Orbits

#### 2.1.1 Introduction to torque

The detailed reaction of a stellar disk to a torque is derived in section 3. Here we introduce a simplified model to establish a framework in which the reaction of the disks to the torques can be studied, and define the terms and symbols we use. Consider a star in the disk on a circular orbit in the $`xy`$-plane around the origin with angular velocity $`\omega (r)=v_c(r)/r`$ (see Figure 1). We apply a torque around the $`x`$-axis by accelerating the star in the $`y`$ and $`z`$ directions around the $`x`$-axis with an antisymmetric pseudo-tensor $`T_{ij}`$. If the acceleration in the $`i`$ direction is $`a_i`$ and the $`j`$-position of the star is $`r_j`$ then the applied acceleration is

$$a_i=\underset{j}{}T_{ij}r_j$$
(1)
where

$$𝑻=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& \tau (r)\\ 0& \tau (r)& 0\end{array}\right).$$

This is an angular acceleration (or equivalently a specific torque) around the $`x`$-axis of magnitude $`a/r=\tau (r)`$. Since the angular acceleration is the time derivative of the angular velocity, the $`x`$-component of the angular velocity increases in response to the torque, and the angular momentum of the star will slew toward the $`x`$ axis, tilting the orbit toward the $`x`$-axis.

#### 2.1.2 Direction of tilt

It is worth taking a minute here to elaborate on the point that the orbit tilts toward the $`x`$ axis. One might think that a torque about the $`x`$ axis would tilt the plane of the orbit around the $`x`$ axis, rather than toward it. Indeed, this is what happens initially, and seems to be the thinking behind attempts to explain the Milky Way warp using the Magellanic Clouds (e.g. Weinberg, 1998; Tsuchiya, 2002). However, García-Ruiz et al (2000) examined the Lagrangian of the system and conducted N-body experiments to demonstrate that a warp caused by the Magellanic Clouds would be perpendicular to the observed Milky Way warp.

This result can be understood by examining the equations of motion of the system for a particle on a circular orbit in the $`xy`$ plane about the $`z`$-axis:

$$\frac{d^2𝒓}{dt^2}=\frac{1}{r}\frac{\varphi }{r}𝒓+𝝉\times 𝒓.$$
(2)
If the untorqued orbit has an angular frequency $`\omega `$, and the torque is directed around the $`x`$-axis, then the equations of motion for the three Cartesian coordinates are

$`{\displaystyle \frac{d^2x}{dt^2}}`$ $`=`$ $`\omega ^2x`$ (3)
$`{\displaystyle \frac{d^2y}{dt^2}}`$ $`=`$ $`\omega ^2y\tau z`$ (4)
$`{\displaystyle \frac{d^2z}{dt^2}}`$ $`=`$ $`\omega ^2z+\tau y.`$ (5)
The $`x`$ motion is decoupled from the other coordinates, and therefore is unaffected by the torque. By substituting equation (4) into equation (5), the system of equations can be converted into a single fourth-order ordinary differential equation for $`z`$ of the form

$$\frac{d^4z}{dt^4}+2\omega ^2\frac{d^2z}{dt^2}+(\omega ^4+\tau ^2)z=0.$$
(6)
Using the ansatz $`z=Ae^{i\mathrm{\Omega }t}`$, the roots of the characteristic equation are

$$\mathrm{\Omega }^2=\omega ^2\pm i\tau .$$
(7)
The two solutions correspond to growing ($``$) and decaying ($`+`$) modes. The growing mode will dominate over time, so we neglect the decaying mode. As $`d^2z/dt^2=\mathrm{\Omega }^2z`$, equations (5) and (7) yield

$$y=iz.$$
(8)
For an oscillating $`y`$, this means that the phase of $`y`$ is a quarter of an orbit ahead of $`z`$, or that the phase of $`z`$ is a quarter of an orbit behind $`y`$. As demonstrated in Figure 2, $`z`$ reaches its maximum a quarter of an orbit after $`y`$, causing the orbit to tilt around the $`y`$-axis toward the $`x`$-axis.

Therefore, we expect that any tilting (and therefore warping, which is a tilt that varies with radius) will occur toward the plane perpendicular to the torque. Since a torque is a transfer of angular momentum, a simple way to understand the result of García-Ruiz et al (2000) is that the galactic disk acquires angular momentum with the same direction as the orbital angular momentum of the satellite, and therefore it will warp toward the plane of the satellite’s orbit, not toward the satellite itself.

#### 2.1.3 Tilting timescale

For an orbit in the $`xy`$-plane, the angular momentum is initially aligned with the $`z`$ axis ($`𝝎=\omega _0(r)`$$`\widehat{𝒛}`$). The torque adds angular momentum in the $`x`$ direction at a constant rate:

$$\frac{d𝝎}{dt}=\tau (r)\widehat{𝒙}.$$

If the radius of the orbit stays constant, then the angular velocity grows as

$$\omega (t)=t\tau (r)\widehat{𝒙}+\omega _0(r)\widehat{𝒛}.$$

As the torque adds angular momentum along the $`x`$ axis, the the orbit tilts toward the $`x`$ axis by an angle $`\theta `$, where

$$\mathrm{tan}\theta =\frac{\omega _x}{\omega _z}=\frac{t\tau (r)}{\omega _0(r)}.$$
(9)
The rate of tilting at small angles is

$$\frac{d\theta }{dt}=\frac{\tau (r)}{\omega _0(r)}$$

with a characteristic timescale to tilt one radian of

$$t_{\mathrm{tilt}}=\frac{\omega _0(r)}{\tau (r)}.$$
(10)
Another motivation for this form is that $`\omega _0`$ is the angular momentum per unit mass and $`\tau `$ is the specific torque, which is the rate of change of the specific angular momentum. Note that even for a flat rotation curve where $`\omega _0=v_c/r`$, $`t_{\mathrm{tilt}}`$ is a monotonically increasing function of $`r`$ if $`\tau (r)`$ falls off faster than $`r^1`$; i.e. for a torque that decreases with radius, the inner regions of the disk should always tilt faster than the outer regions, resulting in a trailing warp.

### 2.2 Expected torques from misaligned flattened halos

There is observational (Sackett, 1999) and theoretical (Katz, 1991; Cole & Lacey, 1996) evidence that dark matter halos are flattened with $`c/a`$ axis ratios ranging as low as 0.5. Cosmological simulations (Frenk et al, 1985; Katz, 1991; Warren et al, 1992; Porciani et al, 2001) show that while the angular momentum and minor axis of a flattened halo will be correlated, there will be significant misalignments between them in a large number of halos. Recent high-resolution simulations suggest that the angular momentum of the baryons may not even be aligned with that of the dark matter (van den Bosch et al, 2002), and that therefore the disks that form from those baryons are usually misaligned with the minor axis of the halo matter distribution.

We calculate the magnitude of the torque for a disk misaligned in a flattened halo with a radial dependence of an NFW profile (Navarro et al, 1996) but flattened along the $`z`$ axis:

$$\rho (x,y,z)=\frac{\rho _0}{\left(m/r_s\right)\left(1+m/r_s\right)^2}$$
(11)
for a modified radius $`m^2=x^2+y^2+z^2/q^2`$ where $`q`$ is the $`c/a`$ axis ratio. The force from this distribution on a given point is calculated using equation (2-88) of Binney & Tremaine (1987). For reference, we also calculate the torques inside flattened isothermal profiles, as in equation (2-54a) of Binney & Tremaine (1987). The torques are calculated by evaluating the forces at opposite points along a fictitious disk centered at the origin and placed at an angle $`\theta `$ to the $`xy`$ plane (see Figure 1). The acceleration at radius $`r`$ due to a torque of magnitude $`\tau `$ is out of the plane, is of opposite sign on opposite sides of the disk, and is of magnitude $`F=\tau r`$. Given the forces on two test points $`F_1`$ and $`F_2`$ at opposite sides of the disk, we take the component of the force out of the plane of the disk $`F_1`$ and $`F_2`$. The torque is

$$\tau =\frac{F_1F_2}{2r}.$$

The angular velocity $`\omega `$ is determined at each point from the radial force: $`\omega =\sqrt{F_r/r}`$ where the radial force $`F_r`$ is defined to be positive when it is directed toward the center of the halo. The average value of $`\omega `$ for the two points across the disk is used. The tilting timescale is given by equation (10).

Figure 3 illustrates the torques derived using this method for halos with a cross-section of properties expected for galactic halos. The fiducial NFW model has flattening $`q=0.7`$, virial velocity $`v_c=175\mathrm{km}\mathrm{s}^1`$ (corresponding to $`v_{\mathrm{rot}}=205\mathrm{km}\mathrm{s}^1`$ at a galactocentric radius of 10 kpc), and a disk at an angle 20° to the symmetry plane of the halo. All NFW halos have concentration parameters $`c_{200}=15`$ (Navarro et al, 1997). The different panels show the effect of varying the disk angle by 10°, the axis ratio from 0.5 to 0.9, and the virial velocity by 75 km s<sup>-1</sup> compared to the fiducial model. The bottom-right panel of Figure 3 shows the torques inside an isothermal profile with a range of axial ratios. The magnitude of the torque in the isothermal case is similar to that in the NFW case, but falls off less rapidly with radius.

Figure 3 demonstrates that the torque from a misaligned flattened halo is significant to its evolution on a cosmological timescale, as the tilting timescale $`t_{\mathrm{tilt}}`$ is less then the Hubble time over the entire length of the disk for all halos modelled. The torques scale approximately as the density, and therefore fall off with radius as $`\tau (r)r^\beta `$ with $`\beta `$ ranging from $`1`$ at small radii to $`3`$ at large radii for NFW halos and equal to $`2`$ for isothermal halos. The typical torque follows the relation

$$\tau (r)=10^{30}\left(\frac{r}{1\mathrm{kpc}}\right)^{2.5}\mathrm{s}^2.$$
(12)
The angular velocity $`\omega (r)r^\alpha `$ where $`\alpha `$ ranges from $`0.5`$ to $`2`$ in NFW halos, and is $`1`$ in isothermal halos. The tilting timescale $`t_{\mathrm{tilt}}`$, as calculated using equation (10), rises almost linearly with radius; the timescale is shorter in the inner regions of the disk, and therefore the inner disk tilts faster than the outer disk in response to torques from misaligned halos. For example, in an isothermal halo $`\tau (r)r^2`$ while $`\omega (r)r^1`$ so $`t_{\mathrm{tilt}}=\omega (r)/\tau (r)r`$. The torque increases more rapidly toward the center than does the disk’s ability to resist the torque due to its angular momentum.

The further the halos and disks are misaligned, the stronger the torque is; however, the torque profile is mostly unchanged for angles beyond 20°. Increasing the virial velocity, and therefore the mass of the halo, increases the magnitude of the torque. Because the concentration $`c_{200}`$ is held constant amongst these models, changing the virial velocity also changes the scale radius $`r_s`$, which can be seen in the shifting of the radius of the knee in the profiles of the top-right panel of Figure 3. The flattening of the halo has a large effect on the magnitude of the torque, with torques strengthening as the halo departs further from spherical symmetry. The torques fall off more slowly in the isothermal halos than in the NFW halos, but are of similar magnitude over most of the disk radius.

## 3 Orbits in Torqued Disks

### 3.1 Reaction of a massless disk

The results of section 2.1 assume that the torque does not change the radius of the orbit. This is a poor approximation, since the torque does work on the star whenever the torque is not orthogonal to the orbital angular momentum. In this section we do a more thorough analysis, taking into account the changing radius of the orbit.

We define $`\theta `$ as in section 2.1, and $`\varphi `$ as the azimuthal angle around the orbit (see Figure 1). Based on the results of section 2.2, the torque can be expressed as a power law:

$$𝝉=\tau _0\left(\frac{r}{r_0}\right)^\beta \widehat{𝒙}$$
(13)
which exerts a force

$$\begin{array}{ccc}\hfill 𝑭& =& m𝝉\times 𝒓\hfill \\ & =& m\tau _0r_0\left(\frac{r}{r_0}\right)^{\beta +1}(\mathrm{cos}\varphi \mathrm{sin}\theta \widehat{𝒚}+\mathrm{sin}\varphi \widehat{𝒛}).\hfill \end{array}$$
(14)
This force torques the orbit and alters the angular momentum $`𝑳`$. If the orbit can be considered to be a circular orbit with circular velocity $`v_c(r)`$ (where the radial dependence is a function of the underlying potential), then

$$𝑳=𝒓\times 𝒑=mrv_c(r)\widehat{𝝎}$$
(15)
where $`\widehat{𝝎}`$ is the unit vector in the direction perpendicular to the orbit. $`𝑳`$ has an additional component of order $`mr^2(d\theta /dt)\widehat{𝒚}`$ because the orbit is tilting, but it does not provide an independent constraint, so we neglect it. The time derivative of $`𝑳`$ is

$$\begin{array}{c}\frac{d𝐋}{dt}=\left[mv_c(r)\mathrm{sin}\theta \frac{dr}{dt}+mr\mathrm{sin}\theta \frac{dv_c}{dr}\frac{dr}{dt}+mrv_c(r)\frac{d\theta }{dt}\mathrm{cos}\theta \right]\widehat{𝒙}\\ +\left[mv_c(r)\mathrm{cos}\theta \frac{dr}{dt}+mr\mathrm{cos}\theta \frac{dv_c}{dr}\frac{dr}{dt}mrv_c(r)\frac{d\theta }{dt}\mathrm{sin}\theta \right]\widehat{𝒛}.\end{array}$$
(16)
The change in angular momentum comes from the torque $`𝑵`$ imparted to the system: $`𝑵\mathbf{=}𝒓\mathbf{\times }𝑭`$. Taking $`𝑭`$ from equation (14) and averaging over one orbit to remove the $`\varphi `$ dependence (assuming that $`dtd\varphi `$ over one orbit, i.e. that $`\omega (r)`$ is slowly varying so the $`\varphi `$-averaged torque is equal to the time-averaged torque) gives

$$\begin{array}{c}𝑵=\frac{1}{2}m\tau _0r_0^2\left(\frac{r}{r_0}\right)^{\beta +2}(1+\mathrm{sin}^2\theta )\widehat{𝒙}\\ +\frac{1}{2}m\tau _0r_0^2\mathrm{sin}\theta \mathrm{cos}\theta \left(\frac{r}{r_0}\right)^{\beta +2}\widehat{𝒛}\end{array}$$
(17)
Setting $`d𝑳/dt=𝑵`$ and solving for $`dr/dt`$ and $`d\theta /dt`$,

$`{\displaystyle \frac{dr}{dt}}`$ $`=`$ $`\tau _0r_0^2\left({\displaystyle \frac{r}{r_0}}\right)^{\beta +2}\left[v_c(r)+r{\displaystyle \frac{dv_c}{dr}}\right]^1\mathrm{sin}\theta `$ (18)
$`{\displaystyle \frac{d\theta }{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\tau _0r_0}{v_c(r)}}\left({\displaystyle \frac{r}{r_0}}\right)^{\beta +1}\mathrm{cos}\theta .`$ (19)
We could not find an analytic solution to this coupled pair of differential equations, so they need to be solved numerically. However, if $`dr/dt`$ is small then equation (19) can be solved at a fixed radius. The solution is

$$\theta (r,t)=2\mathrm{cos}^1\left[\frac{1}{\sqrt{2}}\frac{1+e^{t/t_0(r)}}{\sqrt{1+e^{2t/t_0(r)}}}\right]$$
(20)
where the timescale at a given radius is

$$t_0(r)=\frac{2v_c(r)}{\tau _0r_0}\left(\frac{r}{r_0}\right)^{\beta 1}.$$
(21)
$$=2\frac{\omega _0(r)}{\tau (r)}$$
(22)
and the circular velocity $`v_c(r)`$ depends on the halo potential.

Comparison of equations (22) and (10) reveals that the main effect of allowing a second degree of freedom, by permitting $`r`$ to vary, is to increase the tilting timescale by a factor of two. This factor is also evident from a comparison of the dashed and the solid line in Figure 4. Theoretical studies of warps which artificially restrict this degree of freedom, for example by approximating the disk of stars by solid rings of constant radius, systematically overestimate the effect of the torques on the disk by a factor of two.

### 3.2 Comparison with simulations

The simplified model of equation (9), the detailed derivation (20), and the numerical solution to equations (18) and (19), provide three predictions of how a torqued orbit should evolve. To see how well these predictions work, we compare them to a simulated disk of massless particles subject to a torque that takes all the forces explicitly into account.

We perform a simulation of a disk of 16000 massless tracer particles placed in 800 concentric rings at radii spaced 0.05 kpc apart, each containing 20 particles. They are evolved for 1 Gyr using a predictor-corrector code under the influence of an NFW potential ($`c_{200}=15`$, $`r_{200}=200\mathrm{kpc}`$) and a torquing force of the form (12). We cap the torque at $`\tau =\tau _0`$ inside $`r_0`$ in order to prevent extremely small timesteps for particles near the center where the torquing force diverges. The particles were then binned into spherical shells 3 kpc wide. We calculated the moment of inertia tensor $`I_{ij}=_km_kx_{i,k}x_{j,k}`$ for the particles within each shell and diagonalized it to find the direction of the minor axis (the axis with the smallest moment of inertia). The minor axis was found to be tilted toward the $`x`$-axis by an angle $`\theta `$ that varied from bin to bin.

We also evolved the coupled set of differential equations (18) and (19) forward 1 Gyr with the same values of $`\beta `$ and $`\tau _0(r_0)`$ as above using a fourth-order Runge-Kutta integrator for a set of initial conditions $`\theta =0`$, $`r=0,0.5,1,1.5\mathrm{}100\mathrm{kpc}`$ representing an initially flat disk.

Figure 4 compares the simulation (filled diamonds) with the numerical solution to equations (18) and (19) (open squares), the simplified model (9) (dashed line), and the detailed analytic prediction (20) (solid line). It is apparent that at most radii, (20) is a good approximation to the simulation — much better than (9), demonstrating the importance of including the radial drift of the orbits. At small tilt angles, the difference between the solid line (20) and the dashed line (9) is exactly the factor of two due to the extra degree of freedom derived in Section 3.1. At small radii ($`r5\text{kpc}`$), the solutions diverge: the numerical solution contains rings of the same radius but different tilts while the simulation bins the particles by radius and therefore averages particles in the same shell together; more importantly, the assumption that $`dr/dt`$ is small is violated. However, in the innermost region of a halo the symmetric bulge will dominate the potential, and therefore the torque power law is likely invalid in any case. For these reasons, we ignore the innermost 2 kpc in our analysis and consider equation (20) to be a good description of the disk everywhere else.

## 4 Self-gravitating Disks and Warps

Figure 4 suggests that disks should have continuously curved warps in their inner regions, eventually becoming flat in the outer regions where $`t_{\mathrm{tilt}}`$ is small. Real warped galaxies appear quite the opposite; they are solid out to approximately the Holmberg radius (half of the diameter of the $`I=26.5\mu _{\mathrm{pg}}`$ isophote) and warped beyond that (Briggs, 1990). The most important difference between real galactic disks and those used in section 3 is that the self-gravity of real galactic disks is important; galaxies have mass.

Although the dark matter halo must dominate the potentials of disk galaxies at large radii, there is considerable evidence that the self-gravity of the disk is important to its dynamics over much of its area. Broeils & Courteau (1997) show that “maximal disk” models, where the mass-to-light ratios of the disks $`(M/L)`$ are as large as is consistent with the rotation curves, provide good fits to most rotation curves. Sackett (1997) demonstrates that the Milky Way is also consistent with having a maximal disk. Although Courteau & Rix (1999) have used the residuals of the Tully-Fisher relation (Tully & Fisher, 1977) to estimate that disks provide only 40% of the dynamical mass at 2.2 exponential scale lengths, and de Blok & McGaugh (1997) have demonstrated that the rotation curves of low surface brightness galaxies imply that they are substantially sub-maximal, the prevalent phenomena of spiral arms and bars, which are disk instabilities, require that in many disks, the self-gravity of the disk must be important to its dynamics (Athanassoula et al, 1987).

The self-gravity of the disk acts to keep the disk flat. Ostriker & Binney (1989) examined the effect of a slewing disk potential on a set of self-gravitating rings and found that regions of high surface density react like a solid body. Since the surface density of the disk is highest in the central regions, the central parts of the disk will be kept locally flat, resulting in disks that closer resemble observed warped galaxies. Using numerical N-body simulations, we investigate this effect.

López-Corredoira et al (2002a) take this into account by developing a form for the internal torque $`\tau _{\mathrm{int}}`$. Here we do not use their $`\tau _{\mathrm{int}}`$, but rather integrate the orbits numerically. While this does not provide us with an explicit set of differential equations for $`\theta (r,t)`$, it frees us from such assumptions as the constancy of $`\omega (r)`$ (which we saw was important in Section 3.1), the existence of an equilibrium configuration, and the lack of a dark matter halo.

### 4.1 The simulations

To evaluate how the mass of the disk affects the results of section 3, we performed numerical N-body simulations of massive Milky Way type galactic disks. Equilibrium disk models of disk mass $`1\times 10^{10}\mathrm{M}_{\mathrm{}}`$, $`3\times 10^{10}\mathrm{M}_{\mathrm{}}`$, and $`5.6\times 10^{10}\mathrm{M}_{\mathrm{}}`$, scale length $`r_d=3.5\mathrm{kpc}`$ and vertical scale height $`h_z=325\mathrm{pc}`$ in a static spherically-symmetric NFW halo potential (Navarro et al, 1997) with concentration parameter $`c_{200}=15`$ and virial velocity $`v_{200}=175\mathrm{km}\mathrm{s}^1`$ were constructed using the method of Hernquist (1993) and then allowed to relax under the force of gravity until they appeared to be in equilibrium. The halo masses were varied such that the mass of the disk plus halo was $`4.1\times 10^{11}M_{\mathrm{}}`$ in each case.

The models were evolved using the GRAPESPH code (Steinmetz, 1996), where interparticle gravitational forces are computed using direct summation with GRAPE-3 hardware (Okumura et al, 1993). A plummer softening of 0.3 kpc has been used. The models were evolved for 2 Gyr, which took 5000-7000 timesteps depending on the model. The torque force was applied as an external force as described in section 3.2 with $`\tau _0=10^{30}\mathrm{s}^2`$ at $`r_0=1\mathrm{kpc}`$. Simulations were performed with torque slopes of $`\beta =2.0,`$ $`2.5,`$ and $`3.0`$.

The disks contained 16384 particles. The $`3\times 10^{10}M_{\mathrm{}}`$ simulations were also performed with 32767 particles to see if the resolution was sufficient. The results for the high resolution simulations were identical to those for the lower resolution simulations to well within the angular errors computed in section 4.3 except for in the spherical shells which contained very few particles and had been neglected from the analysis on those grounds.

### 4.2 Warped disks

Figure 5 shows the simulation of a $`3\times 10^{10}M_{\mathrm{}}`$ disk with 32767 particles subject to a $`\beta =2.5`$ torque for 1 Gyr. The main plane of the disk is flat and clearly tilted toward the positive $`x`$-axis, as expected from Section 2.1.2. Beyond 10 kpc, the disk no longer remains flat but warps back toward the original plane. The particles which appear to be filling in the area between the main disk and the warp are projections and are actually in front of or behind the galaxy at large radii.

We binned the simulation particles into spherical shells of thickness 2 kpc and calculated their tilts as in section 3.2. Figure 6 shows a plot of the radial profile of the disk shown in Figure 5, along with the massless simulation and analytic prediction (20) from Figure 4, which correspond to a massless disk in the same torque. There are three important regimes: the central flat disk, the warp, and the effectively massless outer region. The central 11 kpc is all tilted $`17\mathrm{°}`$ from the original plane (we ignore the innermost point because the torque is capped inside $`r_0`$), and corresponds to the visually flat part of the disk. The effect of the disk’s self-gravity is to prevent the inner regions of the disk from tilting as much as they would otherwise, while pulling the outer regions of the disk to larger tilts. At the warp radius $`r_w`$=11 kpc, the disk is no longer flat, and the disk warps back toward the original plane. Finally, from 17 kpc to 23 kpc (the extent of the disk), it follows the analytic prediction, acting like a massless disk. This result agrees qualitatively with that of Ostriker & Binney (1989), who found that regions of high surface density remain flat when torqued, and the position of a warp is determined in part by a drop in surface density.

### 4.3 Warp evolution

We now investigate how the disk evolves over time under the influence of the torque. The general form of the warp is as above, but as the disk evolves under the influence of the torque, the warp radius moves out through the disk at a rate that depends on the mass of the disk. At early times, typically before 200 Myr, the solid region is not well developed and lies below the massless prediction.

We developed an algorithm to automatically detect the warp radius. The particle positions were stored every 40 timesteps. For each of these outputs, the tilt angles were calculated in spherical shells as in Figure 6. We did a bootstrap analysis to estimate the error in these angles: for each spherical shell, we drew 10 random sets from the particles in that shell, calculated the tilt of each bootstrap set, and then used the standard deviation of those angles as an estimate of the error in the tilt of that bin. The error depended on the number of particles in the shell, but was typically less than 1°.

Shells with errors greater than 0.5° were not analyzed, as these generally had few particles and were dominated by numerical noise. For each remaining shell with radius $`r`$ and tilt $`\theta _{\mathrm{sim}}`$, we calculated the deviation $`\mathrm{\Delta }\theta `$ between the simulation tilt $`\theta _{\mathrm{sim}}`$ and the analytic prediction $`\theta _{\mathrm{pred}}(r)`$ from equation (20). The average radius $`r`$ of the radial bin with the highest $`\mathrm{\Delta }\theta `$ was then considered to be the warp radius $`r_w`$. This is the radius at which the difference between the tilt of the massive disk and the tilt of an equivalent massless disk is maximized. For example, in Figure 6, the highest $`\mathrm{\Delta }\theta `$ occurs at $`r=11\text{kpc}`$ where the data point $`\theta _{\text{sim}}=17\mathrm{°}`$ and the predicted line $`\theta _{\text{pred}}=7\mathrm{°}`$. This corresponds exactly to the radial bin at the end of the plateau in Figure 6. Examining a few randomly chosen outputs from each simulation showed that in each case this automated $`r_w`$ agreed with our intuition, except at early times when there was no solid disk and the automated $`r_w`$ was dominated by noise.

Figures 7, 8, and 9 show the evolution of this warp radius for the set of simulations with differing torque slopes $`\beta `$. The squares, circles, and triangles correspond to simulations with disk masses of 1, 3, and $`5.6\times 10^{10}M_{\mathrm{}}`$ respectively. The warp radii appear “quantized” because of the 2 kpc-wide radial bins used to calculate the tilt (see, e.g. Figure 6). One point was taken every $`10^8`$ years, linearly interpolating $`r_w`$ from the neighbouring simulation outputs. The early times before the warp has developed show up as noise in the upper-left region of some of the figures, particularly Figure 9. To quantify the growth of the warp, straight lines were fit through the profiles. While one might think that the lines ought to be constrained to pass through the origin, since there is no warp at $`t=0`$, the data do not support this because the warp instability requires that the surface density of the disk be below a certain critical surface density, as we demonstrate later. We examined Figures 79 by eye and excluded the obviously noisy early times from the fits; in practice, this meant excluding $`t5\times 10^8`$ yr for the $`\beta =3.0`$, $`M_d=1.0\times 10^{10}M_{\mathrm{}}`$ simulation and $`t2\times 10^8`$ yr for the $`\beta =3.0`$, $`M_d=3.0\times 10^{10}M_{\mathrm{}}`$ simulation. We ended the simulations after 2 Gyr, the time when solid region reached the edge of the particle disk in some simulations, at which point the algorithm to find $`r_w`$ just detects the edge of the disk. The fits are shown as the solid, dashed, and dotted lines for the simulations with disk masses of 1, 3, and $`5.6\times 10^{10}M_{\mathrm{}}`$ respectively.

The slope of the linear fit measures how quickly the warp grows. This could depend both on the mass of the disk, which affects the local surface density at the warp radius, and the slope of the torque, which determines the strength of the torque at the warp radius. Figure 10 shows the rate that $`r_w`$ grows for each simulation. There is an apparent trend that the warps in higher mass disks move through the disk faster than in lower mass disks, but the scaling does not appear simple. There is some indication from Figure 10 that the warps from shallower torque laws (which have stronger torques at the warp radii) grow faster than those with steep torques, but the situation is not clear.

The effect of the disk mass on the rate of the warp growth implies that the self-gravity of the disk is important to the formation of the warp. In a massive disk, particles at different radii are coupled to each other gravitationally and act to keep the disk locally flat. This suggests that the crucial factor that determines where the warp develops is the local surface density of the disk. Ostriker & Binney (1989) examined the effect of a slewing disk potential on a set of self-gravitating rings and found qualitatively that regions of high surface density react like a solid body, but that warps can occur where the surface density is lower. Hofner & Sparke (1994) noted that in most cases the group speed of bending waves in a disk of surface density $`\mathrm{\Sigma }(r)`$ and angular rotation velocity $`\omega (r)`$ is

$$c_g=\frac{\pi G\mathrm{\Sigma }(r)}{\omega (r)},$$

and therefore the time for a warp to settle at a given radius would be inversely proportional to the surface density at that radius.

To test this, we translate the warp radii of the simulations into local surface densities using

$$\mathrm{\Sigma }(r_w)=\frac{M_d}{4\pi r_e^2}e^{r_w/r_e}$$
(23)
for a disk of mass $`M_d`$ and exponential scale length $`r_e`$. Figure 11 shows the surface density at the warp radius as a function of time for the three disks of different mass in a $`\beta =2.5`$ torque (i.e. for the same simulations shown in Figure 8), while Figure 12 shows it for a $`3.0\times 10^{10}M_{\mathrm{}}`$ disk in torques with varying $`\beta `$. The surface density at the warp radius falls as the warp moves out through the disk. The local surface densities at the warp *at a given time* are quite similar for all disk masses, and have even smaller scatter for different torque laws. It appears that the local surface density is the important parameter for determining the warp radius at a given time.

The evolution of the warp is well described by a decaying exponential

$$\mathrm{\Sigma }(r_w,t)=\mathrm{\Sigma }_{w_0}e^{t/t_w}$$
(24)
with $`\mathrm{\Sigma }_{w_0}=70M_{\mathrm{}}\mathrm{pc}^2`$ and $`t_w=480`$ Myr, shown as the straight line in Figures 11 and 12. For such a result to hold, the timescale $`t_w`$ should depend only on global properties of the galaxy independent of disk mass or torque; in particular, it can only depend on $`\mathrm{\Sigma }_{w_0}`$, the virial velocity of the galaxy $`v_{200}=175\mathrm{km}\mathrm{s}^1`$, the exponential scale length of the disk $`r_d=3.5\mathrm{kpc}`$, and/or the vertical scale height $`h_z=325\mathrm{pc}`$. An interesting timescale that matches this is the characteristic timescale for bending waves (Hofner & Sparke, 1994) one disk scale length away from $`\mathrm{\Sigma }_{w_0}`$:

$$t_g=\frac{r}{c_g}=\frac{v_{200}}{\pi G\mathrm{\Sigma }_{w_0}e^1}$$
(25)
which is 490 Myr in these cases.

$`\mathrm{\Sigma }_{w_0}`$ is the extrapolation of equation (24) to $`t=0`$, and is the critical surface density above which the disk does not develop a warp. At higher surface densities, the self-gravity of the disk is always sufficient to keep the disk flat. It is interesting that the extrapolation to $`t=0`$, when there is no torque, gives a finite well-defined value for the warp surface density. This is why the linear $`r_w(t)`$ fits do not go through the origin and suggests that at this surface density the disk is marginally unstable to warping. We do not expect to see warps occurring at surface densities higher than $`70M_{\mathrm{}}\mathrm{pc}^2`$. In the Milky Way, the warp begins at or slightly beyond the solar circle (Binney, 1992), while recent Hipparcos determinations of the Milky Way surface density in the solar neighbourhood give $`\mathrm{\Sigma }_0=40M_{\mathrm{}}\mathrm{pc}^2`$ (Crézé et al, 1998). It is encouraging that this is less than $`\mathrm{\Sigma }_{w_0}`$. In external galaxies, warps begin between 25 and 26.5 mag arcsec<sup>-2</sup> (Briggs, 1990), which for a mass to light ratio in $`B`$ of $`1M_{\mathrm{}}/L_{\mathrm{}}`$ gives surface densities of between 1.8–7.0 $`M_{\mathrm{}}\mathrm{pc}^2`$, ignoring projection effects and extinction. This is also consistent with a picture in which warps only can occur at surface densities below $`\mathrm{\Sigma }_{w_0}`$.

## 5 Summary

We studied how a typical galactic disk reacts to torques expected from lying in misaligned dark matter halos. Our main findings are the following:

* Cosmological N-body simulations suggest that galactic disks will be misaligned with the mass distribution of the dark matter halo in which they are embedded. Because of this misalignment, the halo exerts a perpendicular gravitational torque on the disk of the form

$$\tau (r)=\tau _0\left(\frac{r}{1\mathrm{kpc}}\right)^\beta $$

with typical values of $`\tau _0=10^{30}\mathrm{s}^2`$ and $`\beta =2.5`$. This gravitational torque is strong enough to have a significant effect on the entire disk within a Hubble time.
* The timescale for an orbit to tilt in response to this torque rises with radius, so inner portions of the disk will realign themselves first, resulting in a warped disk. For a massless disk of stars in circular orbits, a good estimate of how tilted the disk is at a radius $`r`$ after a time $`t`$ is

$$\theta (r,t)=2\mathrm{cos}^1\left[\frac{1}{\sqrt{2}}\frac{1+e^{t/t_0(r)}}{\sqrt{1+e^{2t/t_0(r)}}}\right]$$

where the timescale $`t_0(r)`$ is defined as

$$t_0(r)=\frac{2v_c(r)}{\tau _0r_0}\left(\frac{r}{r_0}\right)^{\beta 1}$$

and is twice as long as in the case where the orbital radii are not allowed to vary.
* Massive disks depart from this due to the self-gravity of the inner portions of the disk. The disk is kept flatter where the local surface density is high than for a massless disk. The radius inside which the disk is flat (which corresponds to the radius where the warp would be considered to start if the disk were observed) grows with time. The warp grows at a rate of between 3 and 10 kpc Gyr<sup>-1</sup>. The rate depends on the mass of the disk, but is relatively insensitive to the torque parameters.
* More massive disks have faster-growing warps, since only at large radius is the surface density sufficiently low that its self-gravity cannot maintain the flat disk. The surface density at which the warp occurs is well-described by

$$\mathrm{\Sigma }_w(t)=\mathrm{\Sigma }_{w_0}e^{t/t_w}$$

where $`\mathrm{\Sigma }_{w_0}=70M_{\mathrm{}}\mathrm{pc}^2`$ and $`t_w=480`$ Myr, independent of disk mass and torque slope. $`\mathrm{\Sigma }_{w_0}`$ is the maximum surface density at which a warp can occur, and marks the point where a disk is marginally unstable to warping. $`t_w`$ appears to coincide with the timescale of bending waves one disk scale length away from $`\mathrm{\Sigma }_{w_0}`$:

$$t_g=\frac{v_{200}}{\pi G\mathrm{\Sigma }_{w_0}e^1}.$$

In the future, we plan to analyze high resolution cosmological simulations to determine the torques due to realistic dark matter halos. This will enable us to distinguish the torque due to the misaligned disk, that due to substructure in the halo, and that due to the external tidal field by modeling individual halos. A further direction that we have not yet explored is sampling cosmological simulations at regular intervals up to the present day to replace the static torque force used in this work with a torquing field that changes over time in a cosmologically-motivated manner.

Several figures were produced using the WIP package (Morgan, 1995). This work was supported by NSF grants 9807151 and PHY99-0749, and NSERC grant PGSB-233028-2000. MS is an Alfred P. Sloan Fellow and a David and Lucile Packard Fellow.


format: text

text:
# Strategies for High-Throughput, Templated Zeolite Synthesis

## 1 Introduction

Zeolites have found wide industrial applications in processes such as catalysis, molecular sieving, gas separation, and ion exchange. The significant catalytic activity and selectivity of zeolite materials are attributed to the large internal surface area and highly distributed active sites that are accessible through uniformly-sized pores. However, due to the desire to improve the quality of known zeolites that are small or faulted and the continuing quest for new zeolites with different functionalities, novel syntheses of zeolites with various frameworks and chemical compositions have drawn heavy research attention. To date, roughly 150 framework structures have been reported. New structures are added at the pace of roughly 6 per year. The use of organo-cation template molecules to provide structure direction has given rise to breakthroughs in zeolite science and provides much of the impetus for current synthetic efforts . The addition of organic templates such as amines and alkylammonium ions to zeolite synthesis gels can affect the rate at which a particular material is formed and can make new structures or framework chemical compositions accessible. In some cases, a strong correlation between the geometry of the structure directing organic molecule and the zeolite pore architecture is found. The nature and extent of interactions between the organic templates and the inorganic components of a zeolite synthesis gel are key parameters that determine the final zeolite pore architecture .

We are interested in high-throughput synthesis of zeolites. High-throughput, or combinatorial, methods allow for simultaneous creation of a large number of structurally diverse and complex compounds, generalizing the traditional techniques of single compound synthesis. Monte Carlo methods have been proposed and shown to be efficient methods for library design and re-design in both material discovery  and small molecule design . Material discovery deals with continuous variables, such as composition variables and non-composition variables. Small molecule design deals with discrete variables, such as the identities of template and ligand. For templated zeolite synthesis, we have both the continuous zeolite composition and non-composition variables and the discrete identity variables of the component parts of the organic structure-directing molecules. All of these variables affect the final zeolite material in an interrelated way.

We here propose several strategies for templated, high-throughput zeolite synthesis. In section 2 we introduce a new random energy model that serves as a surrogate for experimental measurement of the figure of merit. This random energy model is based upon models used in material discovery and small molecule design. The random energy model is not fundamental to the protocols; it is introduced as a simple way to test, parameterize, and validate the various searching methods. In an experimental implementation, the random energy model will be replaced by the value returned by the screen. In section 3, we introduce several protocols for experimentation. In section 4, we compare results of each protocol and discuss some implications. Finally, we conclude in section 5.

## 2 Random Energy Model

To test proposed protocols in an efficient fashion, we need a model in lieu of the real experimental screening process. That is, we need a model that relates the chemical composition of the zeolite and structure directing agent to the function, or figure of merit, of the material. This model is not essential to the protocols. It is simply a fast and cheap means to test proposed protocols. This model will capture the essence of the physical system and provide validation for the protocols. Random energy models have proven useful for the study of several complex fields, including protein folding , protein molecular evolution , material discovery , and small molecule design .

To consider template-assisted zeolite synthesis, we use a model that combines ideas from both material discovery  and small molecule design . The value of the figure of merit is naturally a sum of a zeolite energy and a well-chosen organic molecular energy:

$$E=E_{\mathrm{zeolite}}+E_{\mathrm{molecule}}$$
(1)
The zeolite framework composition variables $`x_i`$, $`d`$ in number, are certainly key variables in $`E_{\mathrm{zeolite}}`$. The composition variables are specified by the mole fractions $`x_i`$, with $`0x_i1`$, and $`_{i=1}^dx_i=1`$. Typical elemental compositions in zeolite frameworks include silicon, oxygen, aluminum, phosphorus, boron, and germanium, among others. As discussed in the Appendix, the allowed compositions are constrained to lie within a simplex in $`d1`$ dimensions. Several phases will exist for different compositions of the material. The figure of merit will be dramatically different between each of these distinct phases. To mimic this behavior, the composition variables are grouped in the Random Phase Volume Model into phases centered around $`N_x`$ points $`𝐱_\alpha `$ randomly placed within the allowed composition range. The phases form a Voronoi diagram, see Figure 1. The model is defined for any number of composition variables, and the number of phase points is defined by requiring the average spacing between phase points to be $`\xi =0.25`$. To avoid edge effects, additional points are added in a belt of width $`2\xi `$ around the simplex of allowed compositions. The figure of merit should change dramatically between composition phases. Moreover, within each phase $`\alpha `$, the figure of merit should also vary with $`𝐱𝐱_\alpha `$ due to crystallinity effects such as crystallite size, intergrowths, defects, and faulting. Non-composition variables should also affect the measured figure of merit. The non-composition variables are denoted by the $`b`$ variables $`z_i`$. Without loss of generality, each variable ranges $`1z_i1`$. How the figure of merit changes with the non-composition variables should depend on the value of the composition variables. The non-composition variables fall within $`N_z`$ non-composition phases. There are a factor of 10 fewer non-composition phases than composition phases. When $`𝐱`$ is in composition phase $`\alpha `$ and non-composition phase $`\gamma `$, the zeolite energy is given by

$$E_{\mathrm{zeolite}}(𝐱,𝐳)=H(𝐱𝐱^\alpha ,\{G^\alpha \})+\frac{1}{2}H(𝐳,\{G^\gamma \})$$
(2)
where $`\{G^\alpha \}`$ and $`\{G^\gamma \}`$ are phase-dependent random Gaussian variables with zero mean and unit variance. The $`\{G^\alpha \}`$ are different in each composition phase, and the $`\{G^\gamma \}`$ are different in each non-composition phase. The form of the function $`H`$ is

$$H(w_1,\mathrm{},w_n,\{G\})=G_0+\sigma _H\underset{k=1}{\overset{q}{}}\underset{\stackrel{i_1+\mathrm{}+i_n=k}{i_1,\mathrm{},i_n0}}{}f_{i_1\mathrm{}i_n;k}\xi _H^kG_{i_1\mathrm{}i_n}w_1^{i_1}\mathrm{}w_n^{i_n}$$
(3)
The degree of the polynomial is chosen as $`q=6`$. The set of coefficients $`\{G\}`$ includes $`G_0`$ and $`G_{i_1\mathrm{}i_n}`$. The constant symmetry factors $`f_{i_1\mathrm{}i_n;k}`$ are given by

$$f_{i_1\mathrm{}i_n;k}=\frac{k!}{i_1!\mathrm{}i_n!}$$
(4)
The scale factor $`\xi _H`$ is chosen so that each term in the multinomial contributes roughly the same amount in the root-mean-square sense. For the composition variables $`\xi _H=\xi /2`$, and for the non-composition variables $`\xi _H=(z^6/z^2)^{1/4}=(3/7)^{1/4}`$. The $`\sigma _H`$ are chosen so that the multinomial terms contribute $`40\%`$ as much as the corresponding constant, phase terms, $`G_0`$ in eq 3, in the root-mean-square sense.

The organic structure-directing molecules are uniquely characterized by composition, such as the identity of the template and ligands, as shown in Figure 1. Note that *template* here denotes the framework of the structure-directing molecule, onto which ligands are attached, not the entire structure-directing molecule. For simplicity, we consider the small molecule to consist of one template and six binding ligands, each drawn from a fragment library. To quantitatively describe each template and ligand, six weakly correlated descriptors are used. Their numeric values specify the characteristics of each template or ligand.

We again use a random energy model for the molecules to account for the contributions to structure direction arising from interactions between the zeolite and template and from interactions between the zeolite and each of the ligands. In addition, synergistic effects between the ligands and template are incorporated. Consider, for example, a molecule made from template number $`m`$ from the template library and ligands number $`s_1,\mathrm{},s_6`$ from the ligand library. The template is characterized by six descriptors, $`D_1^{(m)},\mathrm{},D_6^{(m)}`$. Similarly, each ligand is characterized by six descriptors, $`d_1^{(l_i)},\mathrm{},d_6^{(l_i)}`$. We denote the template contribution to binding by $`E_\mathrm{T}`$ and the ligand contributions by $`E_\mathrm{L}`$. We denote the contribution due to synergistic ligand-ligand interactions by $`E_{\mathrm{LL}}`$ and the contribution due to synergistic template-ligand interactions by $`E_{\mathrm{TL}}`$. The total molecular contribution to the figure of merit is, then,

$$E_{\mathrm{molecule}}=E_\mathrm{L}+E_\mathrm{T}+E_{\mathrm{LL}}+E_{\mathrm{TL}}$$
(5)
Each of these terms is given in the form of a random polynomial:

$`E_\mathrm{L}(𝐱,𝐳)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{6}{}}}F[d_1^{(l_i)},\mathrm{},d_6^{(l_i)},\{G_\mathrm{L}(𝐱,𝐳)\}],G_\mathrm{L}(𝐱,𝐳)=\lambda _{x,1}G_\mathrm{L}^\alpha +{\displaystyle \frac{1}{2}}\lambda _{z,1}G_\mathrm{L}^\gamma `$ (6)
$`E_\mathrm{T}(𝐱,𝐳)`$ $`=`$ $`F[D_1^{(m)},\mathrm{},D_6^{(m)},\{G_\mathrm{T}(𝐱,𝐳)\}],G_\mathrm{T}(𝐱,𝐳)=\lambda _{x,2}G_\mathrm{T}^\alpha +{\displaystyle \frac{1}{2}}\lambda _{z,2}G_\mathrm{T}^\gamma `$ (7)
$`E_{\mathrm{LL}}(𝐱,𝐳)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{6}{}}}h_i(𝐱,𝐳)F[d_{j_1}^{(l_i)},d_{j_2}^{(l_i)},d_{j_3}^{(l_i)},d_{j_4}^{(l_{i+1})},d_{j_5}^{(l_{i+1})},d_{j_6}^{(l_{i+1})},\{G_{\mathrm{LL}}(𝐱,𝐳)\}],`$ (8)
$`G_{\mathrm{LL}}(𝐱,𝐳)=\lambda _{x,3}G_{\mathrm{LL}}^{i,\alpha }+{\displaystyle \frac{1}{2}}\lambda _{z,3}G_{\mathrm{LL}}^{i,\gamma }`$
$`E_{\mathrm{TL}}(𝐱,𝐳)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{6}{}}}h_i(𝐱,𝐳)F[d_{k_1}^{(l_i)},d_{k_2}^{(l_i)},d_{k_3}^{(l_i)},D_{k_4}^{(m)},D_{k_5}^{(m)},D_{k_6}^{(m)},\{G_{\mathrm{TL}}(𝐱,𝐳)\}],`$ (9)
$`G_{\mathrm{TL}}(𝐱,𝐳)=\lambda _{x,4}G_{\mathrm{TL}}^{i,\alpha }+{\displaystyle \frac{1}{2}}\lambda _{z,4}G_{\mathrm{TL}}^{i,\gamma }`$
where individual descriptors from each template and ligand are given in corresponding libraries. The building block for our random energy model as a function of those descriptors is

$$F(w_1,\mathrm{},w_n,\{G\})=\underset{k=0}{\overset{q}{}}\underset{\stackrel{i_1+\mathrm{}+i_n=k}{i_1,\mathrm{},i_n0}}{}f_{i_1\mathrm{}i_n;k}^{\frac{1}{2}}\xi _F^kG_{i_1\mathrm{}i_n}w_1^{i_1}\mathrm{}w_n^{i_n}$$
(10)
The $`\{G\}`$, again, will be composed of a composition and a non-composition piece, as indicated in eqs 69. We use $`q=6`$ and $`n=6`$. Since the values of the descriptors, $`w_i`$, will be drawn from a Gaussian random distribution of zero mean and unit variance, we set the scaling factor, $`\xi _F`$, by

$$\xi _F=\left(\frac{w^q}{w^2}\right)^{\frac{1}{q2}}=\left(\frac{q!}{(q/2)!2^{q/2}}\right)^{\frac{1}{q2}}$$
(11)
In general the polynomial coefficients, $`\{G\}`$ in eq (6-9), are functions of both zeolite variables $`𝐱`$ and $`𝐳`$. However, for simplicity, we assume that the coefficients are phase-dependent only. Both types of zeolite phases have a significant impact on the energy function. The variables $`\{G_\mathrm{L}^\alpha \},\{G_\mathrm{L}^\gamma \},\{G_\mathrm{T}^\alpha \}`$, $`\{G_\mathrm{T}^\gamma \},\{G_{\mathrm{LL}}^{i,\alpha }\},\{G_{\mathrm{LL}}^{i,\gamma }\},\{G_{\mathrm{TL}}^{i,\alpha }\}`$, and $`\{G_{\mathrm{TL}}^{i,\gamma }\}`$ are fixed random Gaussian variables with zero mean and unit variance. These coefficients are different in each phase. The structural constant $`h_i`$ indicates the strength of the interaction at structure directing site $`i`$. This interaction function depends on the zeolite variables $`𝐱`$ and $`𝐳`$ as

$$h_i(𝐱,𝐳)=H(𝐱𝐱^\alpha ,\{G_i^\alpha \})+\frac{1}{2}H(𝐳,\{G_i^\gamma \})$$
(12)
As before $`\{G_i^\alpha \}`$ and $`\{G_i^\gamma \}`$ are sets of phase-dependent fixed random Gaussian variables with zero mean and unit variance. The parameter $`\sigma _H`$ is adjusted so that the multinomial terms contribute roughly the same as the constant, phase terms on average. In $`E_{\mathrm{LL}}`$, only synergistic interactions between neighboring ligands are considered, and it is understood that $`l_7`$ refers to $`l_1`$ in eq 8. In principle, the polynomial in eq 8 could be a function of all 12 descriptors of both ligands. We assume, however, that important contributions come from interactions among three randomly chosen distinct descriptors of ligand $`l_i,d_{j_1}^{(l_i)},d_{j_2}^{(l_i)},`$ and $`d_{j_3}^{(l_i)},`$ and another three randomly chosen distinct descriptors of ligand $`l_{i+1},d_{j_4}^{(l_{i+1})},d_{j_5}^{(l_{i+1})},`$ and $`d_{j_6}^{(l_{i+1})}`$. Similarly, we assume that template-ligand contributions come from interactions between three randomly chosen distinct descriptors of the ligand, $`d_{k_1}^{(l_i)},d_{k_2}^{(l_i)},`$ and $`d_{k_3}^{(l_i)},`$ and another three randomly chosen distinct descriptors of the template, $`D_{k_4}^{(m)},D_{k_5}^{(m)},`$ and $`D_{k_6}^{(m)}`$. Both $`j_i`$ and $`k_i`$ are descriptor indices ranging from 1 to 6. We assume that these indices depend only on the template, since the molecular details of the zeolite have been suppressed. The constants $`\lambda _{x,i}`$ are adjusted so that the synergistic terms in compositional phases will contribute on average as $`E_\mathrm{L}:E_\mathrm{T}:E_{\mathrm{LL}}:E_{\mathrm{TL}}=1:1:2:1.2`$. Furthermore, we adjust the $`\lambda _{x,i}`$ so that the total molecular contribution from the interaction with the compositional variables accounts for roughly 15% of the constant, compositional phase term in $`E_{\mathrm{zeolite}}`$ in the root-mean-square sense. We similarly adjust the values of $`\lambda _{z,i}`$ in non-compositional phases so that the contribution from the interaction with the non-composition variables is also 15%.

Finally, we minimize the energy returned by the model. That is, we sample the figure of merit, $`E`$, by $`\mathrm{exp}(\beta E)`$.

## 3 Protocols for high-throughput experimentation

As in the small molecule case , we first build the template and ligand libraries. We denote the size of the template library by $`N_\mathrm{T}=15`$ and the size of the ligand library by $`N_\mathrm{L}=1000`$. In a real experiment, the six descriptors would then be calculated for each template and ligand. In the simulated experiment, the values of the six descriptors of each ligand and template are extracted from a Gaussian random distribution with zero mean and unit variance. In the simulated experiment, we also associate two sets of random interaction descriptor indices to each template for the interaction terms in eq 8 and 9.

For comparison, we fix the total number of samples to be synthesized at 100000 for all protocols. That is, we keep roughly the same experimental cost.

We consider several single pass protocols. Current high-throughput experiments uniformly tend to perform a single-pass, grid search on all continuous variables. To mimic current experiments, we use such a grid search on the composition and non-composition variables. For the discrete molecular variables, we instead randomly pick one template and six ligands from the corresponding libraries. In another single-pass protocol, which we call random, we search randomly over all variables, *i.e.* we choose the composition variables and non-composition variables at random as well.

We also consider several multi-pass protocols. Unlike single-pass protocols, multi-pass protocols allow us to learn about the system as the experiment proceeds. In analogy to Monte Carlo methods, each round of combinatorial chemistry experiments corresponds to a move in a Monte Carlo simulation. Instead of tracking one system with many configurational degrees of freedom, however, many samples are tracked. This results in a rather diverse population of samples and increases the opportunity for a few zeolite samples to survive more elaborate tests for application performance, tests that are only roughly correlated with the primary screen. The existence of these secondary and primary screens is the reason why sampling the figure of merit to find several promising compounds is so important. We use a Monte Carlo protocol so as to maintain diversity in the samples. Genetic algorithms, while powerful at local optimization, reduce the diversity of the population in the cross-over step and tend to fail diversity tests . To be consistent with current experiments, we prefer to synthesize 1000 samples simultaneously and perform Monte Carlo moves on every sample for 100 rounds. The initial sample configurations are assigned by the random protocol.

We change the zeolite variables as well as the components of the structure directing agent at each round of the Monte Carlo protocol. In the simplest approach, we perturb both the composition variables and the non-composition variables of the zeolite part around the original value using the traditional Metropolis-type method. For the non-composition variables, $`𝐳`$, periodic boundary conditions are used. For the composition variables, $`𝐱`$, reflecting boundary conditions are used, as discussed in the Appendix. For the molecular part, either the template is changed with probability $`p_{\mathrm{template}}`$, or one of the six ligands is picked randomly to change in each Monte Carlo move. The new configurations are updated according to the acceptance rule at $`\beta `$, the inverse of the protocol temperature. All the samples are sequentially updated in one Monte Carlo round. For the molecular part, either a simple Metropolis method or a biased Monte Carlo method can be incorporated. In the simple Metropolis method, the new template or ligand is selected at random with a uniform probability from the corresponding libraries. The final acceptance rule for an unbiased protocol, which is applied after the zeolite and molecular variables have been modified, is

$$\mathrm{acc}(\mathrm{o}\mathrm{n})=\mathrm{min}[1,\mathrm{exp}(\beta \mathrm{\Delta }E)]$$
(13)
The biased Monte Carlo schemes allow us to generate a trial configuration with a probability that depends on the potential energy of that configuration. These schemes have been proven to be highly successful in small molecule high-throughput experiment design , where biased energy forms are constructed from either theory or pre-experiments. In the present case, since the total molecular energy is a function of both the composition and non-composition variables, it is not feasible to construct the bias from pre-experiments. We, therefore, use the theoretical bias directly from the random energy model. For a zeolite at compositional phase $`\alpha `$ and non-compositional phase $`\gamma `$, the bias energy form for template $`m`$ is

$$E^{(m)}(𝐱,𝐳)=F[D_1^{(m)},\mathrm{},D_6^{(m)},\{G_\mathrm{T}(𝐱,𝐳)\}],G_\mathrm{T}(𝐱,𝐳)=\lambda _{x,2}G_\mathrm{T}^\alpha +\frac{1}{2}\lambda _{z,2}G_\mathrm{T}^\gamma $$
(14)
And the bias for ligand $`i`$ is

$$e^{(i)}(𝐱,𝐳)=F[d_1^{(i)},\mathrm{},d_6^{(i)},\{G_\mathrm{L}(𝐱,𝐳)\}],G_\mathrm{L}(𝐱,𝐳)=\lambda _{x,1}G_\mathrm{L}^\alpha +\frac{1}{2}\lambda _{z,1}G_\mathrm{L}^\gamma $$
(15)
These biases tend to exhibit a large gap between a few dominant templates and ligands and the rest. To ensure the participation of more ligands and templates in the strategy, we introduce cutoff energies for the ligand and template, $`e_\mathrm{c}`$ and $`E_\mathrm{c}`$. We choose $`e_\mathrm{c}`$ to be the 21st lowest ligand energy and $`E_\mathrm{c}`$ to be the 4th lowest template energy. With this cutoff energy, the biased energy, $`e_\mathrm{b}^{(i)}`$, for the $`i^{th}`$ ligand becomes

$$e_\mathrm{b}^{(i)}(𝐱,𝐳)=\{\begin{array}{cc}e^{(i)}(𝐱,𝐳)\hfill & \text{if }e^{(i)}>e_\mathrm{c}\hfill \\ e_\mathrm{c}(𝐱,𝐳)\hfill & \text{otherwise}\hfill \end{array}$$
(16)
And the biased energy, $`E_\mathrm{b}^{(m)}`$, for the $`m^{th}`$ template become

$$E_\mathrm{b}^{(m)}(𝐱,𝐳)=\{\begin{array}{cc}E^{(m)}(𝐱,𝐳)\hfill & \text{if }E^{(m)}>E_\mathrm{c}\hfill \\ E_\mathrm{c}(𝐱,𝐳)\hfill & \text{otherwise}\hfill \end{array}$$
(17)
In one Monte Carlo move, if the $`i^{th}`$ ligand in the molecule is to be changed at new zeolite values $`𝐱^{}`$ and $`𝐳^{}`$, the biased probability for selecting another ligand $`l_i^{}`$ from the ligand library is

$$f[E(\mathrm{n})]=\frac{\mathrm{exp}[\beta e_\mathrm{b}^{(l_i^{})}(𝐱^{},𝐳^{})]}{_{j=1}^{N_\mathrm{L}}\mathrm{exp}[\beta e_\mathrm{b}^{(j)}(𝐱^{},𝐳^{})]}$$
(18)
For reverse move, we have

$$f[E(\mathrm{o})]=\frac{\mathrm{exp}[\beta e_\mathrm{b}^{(l_i)}(𝐱,𝐳)]}{_{j=1}^{N_\mathrm{L}}\mathrm{exp}[\beta e_\mathrm{b}^{(j)}(𝐱,𝐳)]}$$
(19)
Similarly, if the template is desired to be changed at new zeolite values $`𝐱^{}`$ and $`𝐳^{}`$, the biased probability for selecting template $`m^{}`$ is

$$f[E(\mathrm{n})]=\frac{\mathrm{exp}[\beta E_\mathrm{b}^{(m^{})}(𝐱^{},𝐳^{})]}{_{j=1}^{N_\mathrm{T}}\mathrm{exp}[\beta E_\mathrm{b}^{(j)}(𝐱^{},𝐳^{})]}$$
(20)
For reverse move, we have

$$f[E(\mathrm{o})]=\frac{\mathrm{exp}[\beta E_\mathrm{b}^{(m)}(𝐱,𝐳)]}{_{j=1}^{N_\mathrm{T}}\mathrm{exp}[\beta E_\mathrm{b}^{(j)}(𝐱,𝐳)]}$$
(21)
To satisfy the detail balance, the acceptance rule becomes

$$\mathrm{acc}(\mathrm{o}\mathrm{n})=\mathrm{min}\{1,\frac{f[E(\mathrm{o})]}{f[E(\mathrm{n})]}\mathrm{exp}(\beta \mathrm{\Delta }E)\}$$
(22)
Parallel tempering is known to be a powerful tool for searching rugged energy landscapes . Parallel tempering can be combined with both Metropolis Monte Carlo and biased Monte Carlo. In parallel tempering, the samples are divided into $`k`$ groups. The first group of samples is simulated at $`\beta _1`$, the second group is at $`\beta _2`$, and so on, with $`\beta _1<\beta _2<\mathrm{}<\beta _k`$. At the end of each round, samples in group $`i<k`$ are allowed to exchange configurations with samples in group $`i+1`$ with probability $`p_{\mathrm{ex}}`$. The acceptance rule for a parallel tempering exchange move is

$$\mathrm{acc}(\mathrm{o}\mathrm{n})=\mathrm{min}[1,\mathrm{exp}(\mathrm{\Delta }\beta \mathrm{\Delta }E)]$$
(23)
where $`\mathrm{\Delta }\beta =\beta _i\beta _{i+1}`$ and $`\mathrm{\Delta }E`$ is the difference in energy between the sample in group $`i`$ and the sample in group $`i+1`$ before the exchange is made. It is important to notice that this exchange step is experimentally cost-free. Nonetheless, this step can be dramatically effective at facilitating the protocol to escape from local minimum. The number of groups, the number of samples in each group, the value of $`\beta _i`$, and the exchange probability, $`p_{\mathrm{ex}}`$, are experimental parameters to be tuned.

## 4 Results and Discussion

The size of the library is fixed at $`N_\mathrm{T}=15`$ and $`N_\mathrm{L}=1000`$. We adjust $`\lambda _{x,i}`$ and $`\lambda _{z,i}`$ so that the synergistic terms will contribute on average as $`E_\mathrm{L}:E_\mathrm{T}:E_{\mathrm{LL}}:E_{\mathrm{TL}}=1:1:2:1.2`$ in the composition and non-composition phases. We also want a model in which the contribution from zeolite and molecule will be roughly of the same order for the optimal configurations obtained by grid or random protocols. Due to the complicated molecular terms, we find the parameters so that this occurs by trial and error. We first adjust all the $`\sigma _H`$ values, then we adjust the $`\lambda _{x,i}`$ and $`\lambda _{z,i}`$ in eqs 69. As mentioned, we adjust the $`\lambda _{x,i}`$ and $`\lambda _{z,i}`$ so that the total molecular contribution from either type of phase accounts for roughly 15% of the corresponding constant phase term in $`E_{\mathrm{zeolite}}`$.

To locate optimal parameters for the protocols, we perform a few short pre-experiments. Since the size of the template library is relatively small, the optimal value for the probability of changing a template is $`p_{\mathrm{template}}=0.02`$. The maximum random displacements are $`|\mathrm{\Delta }𝐱|=0.1/\sqrt{d1}`$ and $`|\mathrm{\Delta }𝐳|=0.2`$ in the composition space and non-composition space, respectively. For simple Metropolis Monte Carlo, the optimal inverse protocol temperature is $`\beta =50`$ for X3 and X4 and $`\beta =20`$ for X5. For the biased Monte Carlo schemes, the optimal inverse protocol temperature is $`\beta _b=500`$ for X3 and X4 and $`\beta _b=200`$ for X5. For biased Monte Carlo combined with parallel tempering, it is optimal to have the samples divided into three subsets, with 25% of the population at $`\beta _1=\frac{1}{2}\beta _b`$, 50% at $`\beta _2=\beta _b`$, and 25% at $`\beta _3=2\beta _b`$. The switching probability, $`p_{\mathrm{ex}}`$, is 0.1. Determination of these parameter values corresponds experimentally to gaining familiarity with the protocol on a new system.

All the protocols are tested with increasing numbers of composition and non-composition variables. Results are shown in fig 2. From this figure, one can see that all multi-round Monte Carlo protocols are better than the single-pass protocols such as grid and random. The simple Metropolis method can find optimum samples that are twice as superior as those from the grid or random protocols. Most importantly, the biased Monte Carlo is especially efficient. The optimum sample energies from biased Monte Carlo far exceed the optimum sample energies from simple Metropolis with even 1000 rounds. Biased Monte Carlo is able to find such favorable samples mainly by finding much more favorable designs for the structure-directing agent. Parallel tempering is noticeably effective only for the more complicated X4 and X5 systems. For the relatively simple systems, such as X3, it is more important to keep all 1000 of the samples at the optimum temperature.

## 5 Conclusion

We have presented a model for high-throughput zeolite synthesis. As in previous studies, multi-pass Monte Carlo methods work better than do single-pass protocols. Sophisticated biased Monte Carlo schemes are highly efficient and much better than simple Metropolis Monte Carlo. Parallel tempering is the best method for complicated systems with 5 or more framework chemical compositional variables.

## Acknowledgment

It is a pleasure to acknowledge stimulating discussions with Daniel L. Cox that lead this work. This research was supported by the Alfred P. Sloan Foundation through a fellowship to M.W.D., the Chevron Research and Technology Corporation, and the National Science Foundation.

## Appendix: reflecting boundary conditions

The valid region of composition variables form a simplex, since they are non-negative and sum to unity. For those points near the corner of the allowed simplex, a high percentage of proposed moves will go beyond the allowed region. If those moves are simply rejected, the composition point will get stuck in the corner for a long time, and the limited sampling chances will be wasted. We want an algorithm that satisfies detail balance while avoiding getting stuck in the corners.

The scheme we have chosen ensures that all moves from any point in the simplex lead to new valid points in the simplex. The scheme basically involves many reflections over the $`d2`$ dimensional hyper-planes that define the boundaries of the simplex. Figure 3 shows how this works in the case of $`d=3`$ for a single reflection. We first define the unit normals, $`𝐧_i`$, to each of the faces of the simplex. We also find the constants $`c_i`$ that define the faces by the equation $`𝐱𝐧_i=c_i`$. In the case of figure 3, there are three such unit normals and three such constants. We initially define $`𝐲=𝐱(\mathrm{o})`$. The algorithm consists of the following steps:

1. Test whether the proposed new point, $`𝐱`$, is in the allowed region. If so, define the new point $`𝐱(\mathrm{n})`$ to be $`𝐱`$ and stop.
2. If not, find the face $`i`$ for which the quantity $`t_i=(c_i𝐲𝐧_i)/[(𝐱𝐲)𝐧_i]`$ is minimal, taking into account only those faces for which both $`t_i`$ and $`(𝐱𝐲)𝐧_i`$ are positive.
3. Define $`𝐲𝐲+t_i(𝐱𝐲)`$.
4. Reflect the point through this face by the equation $`𝐱𝐱2(𝐱𝐧_ic_i)𝐧_i`$.
5. Continue with step 1.

This algorithm converges, because it decreases the magnitude of $`𝐱`$ at each step, until $`𝐱`$ is within the allowed region. The algorithm is also reversible. That is, for each forward move $`𝐱(\mathrm{o})𝐱`$ that maps to the new position $`𝐱(\mathrm{n})`$, there is always a reverse move $`𝐱(\mathrm{n})𝐱^{}`$ that maps to the original position $`𝐱(\mathrm{o})`$. Furthermore, when the original move is in the allowed move sphere, $`|𝐱(\mathrm{o})𝐱|<\mathrm{\Delta }x`$, the reverse move is as well, $`|𝐱(\mathrm{n})𝐱^{}|<\mathrm{\Delta }x`$. Given that we impose the detailed balance condition, then, these reflecting boundary conditions provide a valid Monte Carlo move from an old point in the allowed simplex to a new point in the allowed simplex.

These reflecting boundary conditions can be interpreted in a geometrical way. Essentially, the reflections perform a type of billiards in the $`d1`$ dimensional simplex. In each attempted move, we chose a displacement $`𝐱𝐱(\mathrm{o})`$. We can interpret the composition variable as a small ball that moves at a constant speed along this trajectory. The motion is continued until either a boundary is encountered or the entire length of the trajectory has been traveled. If a boundary of the simplex is encountered, the ball reflects off the hyper-plane by Newtonian mechanics, and the trajectory continues along the new, reflected direction. The motion, including possible additional reflections, is continued until the ball has traveled a distance equal to the chosen length of the displacement, $`|𝐱𝐱(\mathrm{o})|`$. The location of the ball at the end of the trajectory then gives the new point $`𝐱(\mathrm{n})`$. The reverse move is given by exactly the reverse motion, which has the same total displacement.