text:
# Surface Outflow in Optically Thick Dust Disks by the Radiation Pressure 1footnote 11footnote 1To appear in Aug. 10 issue of Astrophys. J.

## 1 Introduction

Circumstellar disks around young stars are composed of gas and dust particles, and in general, the dust and gas accretion velocities are not equal. The gas suffers turbulent viscosity and accretes onto the central star on a viscous time-scale of order $`10^610^7`$ yr (Hartmann et al. 1998). The orbital evolution of dust particles, on the other hand, is controlled by gas drag (Adachi, Hayashi, & Nakazawa 1976; Weidenschilling 1977).

Several factors affect the radial motion of particles. First, gas drag in the radial direction entrains particles in the accretion flow. When the gas density is high or when the particle size is small, the particles follow the gas motion, and accrete onto the star on the viscous time-scale. Second, angular momentum is transfered between particles and the gas. The radial gas pressure gradient shifts the angular velocity of the gas slightly from the Keplerian value. Moreover, if starlight can irradiate particles, then radiation pressure causes non-Keplerian deviations in the particle velocity as well. Dust particles can thus gain or lose angular momentum through gas drag and thus move outward or inward, depending on whether they experience a net head-wind or a net tail-wind.

In a previous paper (Takeuchi & Lin 2002, hereafter Paper I), we studied particle motion in optically thick disks, using the assumption that radiation pressure on the dust is negligible. We found that dust particles smaller than $`10\mathrm{\mu m}`$ accrete with smaller velocity than the gas. As the gas accretion proceeds, small dust particles are left behind in the disk and the dust-to-gas ratio gradually increases. Because circumstellar disks associated with young stars are generally optically thick, the majority of dust particles are not directly exposed to the stellar radiation, and the approximation used in Paper I is applicable. There are always, however, some particles near the disk surfaces that are directly irradiated, and the motion of such particles are affected by radiation pressure. Klahr & Lin (2001) and Takeuchi & Artymowicz (2001) studied the motion of particles in optically thin gas disks around young Vega-type stars, where the dust experiences both gas drag and radiation pressure. Because the radiation pressure reduces the angular velocity of the dust below that of the gas, the dust tends to move outward. Thus, in a more refined description of an optically thick disk, one expects that particles in the irradiated surface layer move outward, while beneath the surface layer particles move inward. If the surface outflow dominates the midplane inflow, then there is a net radial expansion that may induce the formation of rings or gaps.

In this paper, we study the surface dust outflow in optically thick disks due to radiation pressure. The mass flux of the surface outflow is calculated, which is compared to that of the inflow inside the disk. We find that the surface outflow is negligible in a minimum-mass-solar-nebula model, but that the surface outflow can dominate inflow at some disk radii in some cases, including situations involving high stellar luminosity or low gas mass. In such disks, the appearance of rings or gaps is expected.

## 2 The Distribution of the Gas and the Dust

### 2.1 The Gas Disk

We consider circumstellar disks around young pre-main-sequence stars. The disks are turbulent and are accreting onto the central stars, but they are assumed to be passive, that is, the disks are primarily heated by stellar radiation only. When turbulent fluctuations are averaged out, the mean gas motion is nearly in hydrostatic equilibrium and the radial accretion flow is much slower than the azimuthal Kepler rotation. In cylindrical coordinates ($`r,\theta ,z`$), the balance of forces in the $`r`$\- and $`z`$-directions are

$$r\mathrm{\Omega }_g^2\frac{GMr}{(r^2+z^2)^{3/2}}\frac{1}{\rho _g}\frac{P_g}{r}=0,$$
(1)
$$\frac{GMz}{(r^2+z^2)^{3/2}}\frac{1}{\rho _g}\frac{P_g}{z}=0,$$
(2)
where $`\mathrm{\Omega }_g`$ is the gas angular velocity, $`M`$ is the central star’s mass, $`\rho _g`$ is the gas density, and $`P_g`$ is the gas pressure. For simplicity, we assume that the gas disk has power-law density and temperature profiles in the radial direction, and is isothermal in the vertical direction. The temperature distribution is written as a power-law of index $`q`$,

$$T(r,z)=T_0r_{\mathrm{AU}}^q,$$
(3)
where the subscript “0” denotes quantities at 1 AU, and a non-dimensional quantity $`r_{\mathrm{AU}}`$ is the radius in AU. The isothermal sound speed is $`c=c_0r_{\mathrm{AU}}^{q/2}`$ and the gas disk scale height is thus defined as

$$h_g(r)\frac{c}{\mathrm{\Omega }_{\mathrm{K},\mathrm{mid}}}=h_0r_{\mathrm{AU}}^{(q+3)/2},$$
(4)
where $`\mathrm{\Omega }_\mathrm{K}=[GM/(r^2+z^2)^{3/2}]^{1/2}`$ is the Keplerian angular velocity and the subscript “mid” indicates the midplane value. From the balance between the stellar gravity and the gas pressure gradient in the $`z`$-direction (eq. ), we have the density distribution

$$\rho _g(r,z)=\rho _0r_{\mathrm{AU}}^p\mathrm{exp}\left(\frac{z^2}{2h_g^2}\right),$$
(5)
where a power-law in $`r`$ is assumed. The surface density of the gas disk is given by

$$\mathrm{\Sigma }_g(r)=_{\mathrm{}}^+\mathrm{}\rho _g𝑑z=\sqrt{2\pi }\rho _0h_0r_{\mathrm{AU}}^{p_s},$$
(6)
where $`p_s=p+(q+3)/2`$. The angular velocity $`\mathrm{\Omega }_g`$ of the gas is slightly different from the Keplerian angular velocity $`\mathrm{\Omega }_\mathrm{K}`$ because of the radial pressure gradient in the gas. From equation (1), we have

$$\mathrm{\Omega }_g^2(r,z)=\left(1\eta \right)\mathrm{\Omega }_\mathrm{K}^2(r,z),$$
(7)
where $`\eta =(r\mathrm{\Omega }_\mathrm{K}^2\rho _g)^1P_g/r`$ is the ratio of the gas pressure gradient to the stellar gravity in the radial direction. As derived in equation (17) in Paper I, $`\eta `$ is written to order $`(h_g/r)^2`$ as

$$\eta =\left(\frac{h_g}{r}\right)^2\left(p+q+\frac{q+3}{2}\frac{z^2}{h_g^2}\right).$$
(8)
The disk has a turbulent viscosity $`\nu `$. Turbulent disk viscosity has been ascribed to various mechanisms, including the convective instability (Lin & Papaloizou 1980) and the magneto-rotational instability (Balbus & Hawley 1991). In this paper, we simply model the viscous effect of turbulence using the so-called $`\alpha `$ prescription (Shakura & Sunyaev 1973),

$$\nu =\alpha ch_g=\alpha c_0h_0r_{\mathrm{AU}}^{q+3/2}.$$
(9)
We adopt the following fiducial parameters: $`M=1M_{\mathrm{}}`$, $`\rho _0=2.83\times 10^{10}\mathrm{g}\mathrm{cm}^3`$, $`h_0=3.33\times 10^2`$ AU, $`p=2.25`$, $`q=0.5`$, and $`\alpha =10^3`$. These values correspond to a disk having a gas mass of $`2.5\times 10^2M_{\mathrm{}}`$ inside 100 AU. The power law index of the surface density is $`p_s=p+(q+3)/2=1.0`$. The temperature distribution is $`T=278r_{\mathrm{AU}}^{1/2}\mathrm{K}`$.

### 2.2 The Dust Disk

#### 2.2.1 Radial Distribution

The radial surface density distribution of the dust is assumed to be proportional to that of the gas disk, so that $`\mathrm{\Sigma }_{d,\mathrm{all}}(r)=f_{\mathrm{dust}}\mathrm{\Sigma }_g(r)`$, where the constant$`f_{\mathrm{dust}}`$ is the dust-to-gas mass ratio. Throughout the paper, we adopt the value $`f_{\mathrm{dust}}=0.01`$ at all radii. The subscript “all” represents contributions from the entire distribution of particle sizes. In this paper, we consider cases where all particles have the same size and also power-law size distributions. For power-law size distributions, the size range is bracketed at the large and small ends, $`[s_{\mathrm{min}},s_{\mathrm{max}}]`$. We adopt a power-law index $`3.5`$, appropriate to the distribution expected following collisional evolution (Hellyer 1970). The surface number density of particles in a size range $`[s,s+ds]`$ at a distance $`r`$ from the star is given by

$$N_\mathrm{\Sigma }(r,s)ds=N_1(r)s_{\mathrm{cm}}^{3.5}ds,$$
(10)
where $`s_{\mathrm{cm}}`$ is the particle size in centimeters. The radial distribution of particles $`N_1`$ is

$$N_1(r)=\frac{3}{8\pi \rho _p}f_{\mathrm{dust}}\mathrm{\Sigma }_g(r)\left(s_{\mathrm{max},\mathrm{cm}}^{1/2}s_{\mathrm{min},\mathrm{cm}}^{1/2}\right)^1(\mathrm{cm})^4,$$
(11)
where $`\rho _p`$ is the physical density of a particle. We adopt $`\rho _p=1.25\mathrm{g}\mathrm{cm}^3`$, $`f_{\mathrm{dust}}=0.01`$, $`s_{\mathrm{min}}=0.01\mathrm{\mu m}`$, and $`s_{\mathrm{max}}=1\mathrm{\mu m}`$.

#### 2.2.2 Vertical Distribution

Dust particles feel the $`z`$-component of the stellar gravity, and thus sediment toward the disk midplane. Small particles ($`s1`$ cm) immediately (within one orbit) reach the terminal velocity $`v_{z,d}`$, where the gravity and the gas drag are in balance. For small particles considered here ($`s1`$ cm), the Epstein’s drag law is applied, because the mean free path of the gas molecules (7 cm at 1 AU) is larger than the particle size. The gas drag force is

$$𝑭_g=\frac{4}{3}\pi \rho _gs^2v_t(𝒗_d𝒗_g),$$
(12)
where $`v_t`$ is the mean thermal velocity of the gas molecules, and $`𝒗_d`$ and $`𝒗_g`$ are the velocities of the dust and the gas, respectively. The stopping time due to the gas drag is $`t_s=m_d|𝒗_d𝒗_g|/|𝑭_g|`$, where $`m_d=\frac{4}{3}\pi \rho _ps^3`$ is the particle mass. We define the non-dimensional stopping time normalized by the orbital time as,

$$T_s=t_s\mathrm{\Omega }_{\mathrm{K},\mathrm{mid}}=\frac{\rho _psv_{\mathrm{K},\mathrm{mid}}}{\rho _grv_t},$$
(13)
where $`v_{\mathrm{K},\mathrm{mid}}`$ is the Keplerian velocity at the midplane. When terminal velocity is achieved, the gas drag balances the $`z`$-component of the stellar gravity, $`m_d\mathrm{\Omega }_\mathrm{K}^2z`$. Using the approximation $`\mathrm{\Omega }_\mathrm{K}\mathrm{\Omega }_{\mathrm{K},\mathrm{mid}}`$, the terminal velocity $`v_{z,d}`$ is

$$v_{z,d}=\mathrm{\Omega }_{\mathrm{K},\mathrm{mid}}T_sz.$$
(14)
The mass flux of particles in the vertical direction is

$$f_{z,\mathrm{sed}}=\rho _dv_{z,d},$$
(15)
where $`\rho _d`$ is the dust mass density. If the gas disk were laminar, the particles would concentrate at the midplane in a time that is short compared to the evolution time of the gas disk (the sedimentation time-scale is about $`10^6`$ yr for $`1\mathrm{\mu m}`$ particles at 10 AU). When the gas disk is turbulent, however, the gas stirs and disperses the particles and prevents sedimentation. The mass flux of dust arising from this turbulent diffusion is assumed proportional to the gradient of the dust concentration $`\rho _d/\rho _g`$ in analogy with molecular diffusion (e.g, Monin & Yaglom, 1971, p. 579; Morfill, 1985, §3.5.1; Cuzzi, Dobrovolskis, & Champney 1993). In the vertical direction, the mass flux $`f_{z,\mathrm{dif}}`$ is

$$f_{z,\mathrm{dif}}=\frac{\rho _g\nu }{\mathrm{Sc}}\frac{}{z}\left(\frac{\rho _d}{\rho _g}\right),$$
(16)
where the Schmidt number, $`\mathrm{Sc}`$, represents the coupling strength between the particles and the gas. For small particles that are well coupled to the gas, $`\mathrm{Sc}`$ approaches unity. For large particles that are decoupled from the gas, Sc becomes infinite. In our standard model, we use $`\mathrm{Sc}=1`$.

In equilibrium, the sedimentation mass flux balances the diffusive mass flux:

$$f_{z,\mathrm{sed}}+f_{z,\mathrm{dif}}=0.$$
(17)
The equilibrium dust density distribution is given by equation (31) in Paper I. For particles in a size range $`[s,s+ds]`$, the density distribution is

$`\rho _d(r,z,s)ds`$ $`=`$ $`\rho _d(r,0,s)\mathrm{exp}[{\displaystyle \frac{z^2}{2h_g^2}}`$
$``$ $`{\displaystyle \frac{\mathrm{Sc}T_{s,\mathrm{mid}}}{\alpha }}(\mathrm{exp}{\displaystyle \frac{z^2}{2h_g^2}}1)]ds,`$
(18)
where $`T_{s,\mathrm{mid}}`$ is the midplane value of the non-dimensional stopping time. The total dust density is

$$\rho _{d,\mathrm{all}}(r,z)=_{s_{\mathrm{min}}}^{s_{\mathrm{max}}}\rho _d(r,z,s)𝑑s.$$
(19)
As discussed in Paper I, the sedimentation of particles is more effective for larger particles and at larger distances from the star. As increasing the distance, the dust disk becomes relatively thinner compared to the gas disk.

### 2.3 Optical Depth from the Star

The disks considered in this paper are optically thick at visible wavelengths even in the vertical direction ($`\tau _{\mathrm{ver}}100`$ at 100 AU). Stellar light cannot penetrate the dust disk. Only particles residing in the thin surface layer are exposed to stellar radiation pressure.

The optical depth from the star to the position $`(r,z)`$ in the dust disk is

$$\tau (r,z)=_{r_{\mathrm{in}}}^r\frac{\kappa }{f_{\mathrm{dust}}}\rho _{d,\mathrm{all}}(r^{},z^{})\sqrt{1+\frac{z^2}{r^2}}𝑑r^{},$$
(20)
where $`z^{}=(z/r)r^{}`$, $`r_{\mathrm{in}}`$ is inner edge of the dust disk, and $`\kappa `$ is the dust opacity in the visible wavelength for an unit mass of the gas. The opacity for a unit mass of the dust is $`\kappa /f_{\mathrm{dust}}`$. In our calculations, we adopt $`r_{\mathrm{in}}=0.1\mathrm{AU}`$. This choice of inner boundary does not affect the integral (20), because the dust disk is flared and the main contribution to the optical depth comes from the dust at relatively large distances from the star. We adopt a dust mass opacity $`\kappa =100\mathrm{cm}^2\mathrm{g}^1`$ for visible wavelengths. For particles smaller than $`100\mathrm{\mu m}`$, the opacity is independent of (the maximum) particle size (see Figs. 5 and 11 in Miyake & Nakagawa 1993). Thus, we need not consider variations in $`\kappa `$ with changes in the particle size.

Figure 1 shows the disk surfaces at which the optical depth to the starlight is unity. Above these illuminated surfaces, the dust particles are directly exposed to the stellar radiation. Figure 1$`a`$ shows the illuminated surfaces of the standard disk models which are composed of single-sized particles ($`0.1`$ to $`10\mathrm{\mu m}`$). As the particle size increases, the level of the illuminated surfaces becomes lower. This size dependence is because larger particles sediment more readily, and thus concentrate near the midplane. The illuminated surfaces have a slightly flatter geometry than the gas disk, because the dust particles at larger distances sediment more effectively. Figure 1$`b`$ shows the illuminated surfaces of the models containing power-law size distributions. The illuminated surfaces lie at roughly three disk scale heights, which are similar to those of the single-sized $`0.1\mathrm{\mu m}`$ particles. The location of the illuminated surface does not change significantly with the maximum size of the particles, because high altitude above the midplane is mostly populated by the smallest particles.

## 3 Outflow in the Surface Layer

Dust particles above the illuminated surface suffer radiation pressure from the central star. The radiation pressure balances a part of the stellar gravity, so that the small particles orbit slower than the gas. Because the gas drag transfers angular momentum from the gas to the particles, the particles move outward, leading to an outflow in the surface layers of dust disks. In this section, we estimate the mass fluxes carried by these outflows.

### 3.1 Outflow Mass Flux

The radiation pressure is expressed through the ratio $`\beta `$ to the stellar gravity. The radial component is

$$F_{\mathrm{rad}}=\beta m_dr\mathrm{\Omega }_\mathrm{K}^2.$$
(21)
The value of $`\beta `$, which is calculated with the Mie’s scattering theory, depends on the particle composition, porosity, shape, and size. In this paper, we adopt $`\beta `$ for “young cometary particles,” which is calculated by Wilck & Mann (1996) for typical interplanetary particles in the solar system and shown in Figure 2. In general, stellar light causes Poynting-Robertson drag on the particles in addition to the gas drag induced by radiation pressure. However, in the gas disks considered here, the gas density is sufficiently large so that gas drag dominates Poynting-Robertson drag (see the discussion in Takeuchi & Artymowicz 2001, §4.2). Poynting-Robertson drag can therefore be neglected.

The equations of motion of a particle are

$$\frac{dv_{r,d}}{dt}=\frac{v_{\theta ,d}^2}{r}(1\beta )\mathrm{\Omega }_\mathrm{K}^2r\frac{\mathrm{\Omega }_{\mathrm{K},\mathrm{mid}}}{T_s}(v_{r,d}v_{r,g}),$$
(22)
$$\frac{d}{dt}\left(rv_{\theta ,d}\right)=\frac{v_{\mathrm{K},\mathrm{mid}}}{T_s}(v_{\theta ,d}v_{\theta ,g}),$$
(23)
where $`v_r`$ and $`v_\theta `$ are the $`r`$\- and $`\theta `$\- components of the velocity, and the subscripts “$`g`$” and “$`d`$” distinguish gas and dust. As in Paper I, we assume $`v_{\theta ,g}v_{\theta ,d}v_{\mathrm{K},\mathrm{mid}}`$ and $`d(rv_{\theta ,d})/dtv_{r,d}d(rv_{\mathrm{K},\mathrm{mid}})/dr=v_{r,d}v_{\mathrm{K},\mathrm{mid}}/2`$. From equation (23) we then have

$$v_{\theta ,d}v_{\theta ,g}=\frac{1}{2}T_sv_{r,d}.$$
(24)
In equation (22), $`dv_{r,d}/dt=O(v_{r,d}^2/r)`$ is neglected if $`(v_{r,d}/v_{\mathrm{K},\mathrm{mid}})^21`$. We also neglect second order terms in $`\eta `$ and $`(v_{\theta ,d}v_{\theta ,g})/v_{\mathrm{K},\mathrm{mid}}`$. To this level of approximation, $`v_{\theta ,d}^2v_{\theta ,g}^2+2v_{\theta ,g}(v_{\theta ,d}v_{\theta ,g})v_\mathrm{K}^2(1\eta )+2v_{\mathrm{K},\mathrm{mid}}(v_{\theta ,d}v_{\theta ,g})`$, and equation (22) becomes

$$(\beta \eta )v_\mathrm{K}^2+2v_{\mathrm{K},\mathrm{mid}}(v_{\theta ,d}v_{\theta ,g})\frac{v_{\mathrm{K},\mathrm{mid}}}{T_s}(v_{r,d}v_{r,g})=0.$$
(25)
Using equation (24) and $`v_\mathrm{K}v_{\mathrm{K},\mathrm{mid}}`$, the particle radial velocity is

$$v_{r,d}(r,z,s)=\frac{T_s^1v_{r,g}+(\beta \eta )v_{\mathrm{K},\mathrm{mid}}}{T_s+T_s^1}.$$
(26)
The value of $`\eta `$ is of order $`(h_g/r)^2`$ and it is smaller than $`0.05`$. Since particles smaller than $`10\mathrm{\mu m}`$ have $`\beta \eta `$ (see Fig. 2), $`\eta `$ in equation (26) can be neglected. Particles larger than $`100\mathrm{\mu m}`$ may have $`\beta `$ smaller than $`\eta `$. For such large particles, the radiation pressure is unimportant, even if they are above the illuminated surface. In this section, we consider only small particles with radii less than $`10\mathrm{\mu m}`$. Such particles have non-dimensional stopping times $`T_s1`$ (see Fig. 2 in Paper I). The radial velocity then reduces to

$$v_{r,d}=v_{r,g}+\beta T_sv_{\mathrm{K},\mathrm{mid}}.$$
(27)
The mass flux of the outflow in the surface layer is

$$F_{\mathrm{sur}}(r)=_{s_{\mathrm{min}}}^{s_{\mathrm{max}}}2_{z_{\mathrm{sur}}}^{\mathrm{}}2\pi r\rho _d(r,z,s)v_{r,d}(r,z,s)𝑑z𝑑s,$$
(28)
where $`z_{\mathrm{sur}}`$ is the height of the illuminated surface. The factor 2 arises from considering both the top and bottom surfaces. In the integration in $`z`$ above the illuminated surface, $`v_{r,d}`$ varies less rapidly than $`\rho _d`$, and we may use its value at the illuminated surface. The mass flux is then expressed as

$$F_{\mathrm{sur}}(r)=2\pi r_{s_{\mathrm{min}}}^{s_{\mathrm{max}}}v_{r,d}(r,z_{\mathrm{sur}},s)\mathrm{\Sigma }_{d,\mathrm{sur}}(r,s)𝑑s,$$
(29)
where $`\mathrm{\Sigma }_{d,\mathrm{sur}}ds`$ is the surface density of dust particles above the illuminated surface in the size range $`[s,s+ds]`$,

$$\mathrm{\Sigma }_{d,\mathrm{sur}}(r,s)=2_{z_{\mathrm{sur}}}^{\mathrm{}}\rho _d(r,z,s)𝑑z.$$
(30)
The surface outflow can clear out the dust within a distance $`r`$ from the star on a time-scale $`t_{\mathrm{out}}M_{\mathrm{dust}}(r)/F_{\mathrm{sur}}(r)`$, where $`M_{\mathrm{dust}}(r)`$ is the dust mass within $`r`$. In our models, the disk mass is dominated by the outermost part of the disk and $`M_{\mathrm{dust}}(r)\pi r^2\mathrm{\Sigma }_{d,\mathrm{all}}(r)`$. The dust clearing (or the disk evolution) time-scale is thus

$$t_{\mathrm{out}}(r)\frac{\pi r^2\mathrm{\Sigma }_{d,\mathrm{all}}(r)}{F_{\mathrm{sur}}(r)}.$$
(31)
While at the surface layer particles flow outward, the particles beneath the surface layer are shielded from stellar irradiation. At such locations, the gas drag causes a net inward flow of particles as shown in Paper I. From equation (26) with $`\beta =0`$, the radial velocity of particles beneath the surface layer is

$$v_{r,d}(r,z,s)=\frac{T_s^1v_{r,g}\eta v_{\mathrm{K},\mathrm{mid}}}{T_s+T_s^1}.$$
(32)
The mass flux is

$$F_{\mathrm{in}}=_{s_{\mathrm{min}}}^{s_{\mathrm{max}}}_{z_{\mathrm{sur}}}^{z_{\mathrm{sur}}}2\pi r\rho _d(r,z,s)v_{r,d}(r,z,s)𝑑z𝑑s.$$
(33)
The time required for the inflow to transport a mass equal to the dust mass within $`r`$ is

$$t_{\mathrm{in}}(r)\frac{\pi r^2\mathrm{\Sigma }_{d,\mathrm{all}}(r)}{|F_{\mathrm{in}}(r)|}.$$
(34)
We define this as the disk evolution time through the inflow. This definition of $`t_{\mathrm{in}}`$ does not mean the time for carrying the total dust mass outside $`r`$ to the inner part, but it is convenient for comparing with $`t_{\mathrm{out}}`$. Comparison of $`t_{\mathrm{out}}`$ and $`t_{\mathrm{in}}`$ simply means comparison of the outflow and inflow fluxes. The disk clearing through the surface outflow occurs effectively only if the surface flux is stronger than the inflow flux, i.e., $`t_{\mathrm{out}}<t_{\mathrm{in}}`$. The comparison of these two time-scales is discussed in §3.3 and §3.4 below.

### 3.2 Evolution Time of the Disks through the Outflow

We first discuss the disk clearing time, $`t_{\mathrm{out}}`$, through the surface outflow. In this subsection, we adopt the assumption that inflow flux near the midplane is negligible and that only the surface flow affects the dust structure within the disk. In §3.3 and §3.4, we relax this assumption and consider more realistic scenarios including the effects of inflow near the midplane.

The particles in the surface layer have outflow velocity, $`v_{r,d}(r,z_{\mathrm{sur}},s)`$. These particles are blown from their orbits to larger radii on a time-scale $`t_b=r/v_{r,d}(r,z_{\mathrm{sur}},s)`$. For disks composed of single-sized particles, $`t_b`$ is plotted in Figure 3. The blow-off time for $`0.1`$ to $`10\mathrm{\mu m}`$ particles is about $`10^2`$ yr at a few AU from the star, and is about $`10^5`$ yr at 100 AU. The blow-off time is much shorter than the inward migration time inside the disk, which is discussed in §3.3 below, because $`\beta \eta `$ and because the stopping time $`T_s`$ at the illuminated surface is much longer than at the midplane.

In our standard model, the dust disk is highly optically thick. Even in the vertical direction at 100AU, the optical depth is $`\tau _{\mathrm{ver}}100`$, so the layer above the illuminated surface is very thin. Figure 4 shows the ratio of the surface density of the layer to that of the whole disk ($`\mathrm{\Sigma }_{d,\mathrm{sur},\mathrm{all}}/\mathrm{\Sigma }_{d,\mathrm{all}}`$). As the distance from the star increases, the mass ratio increases and the surface layer becomes thicker. However, even at 100 AU, the mass of the surface layer is less than $`10^3`$ times the dust disk mass. In the case of single-sized particles (Fig. 4$`a`$) the mass of the surface layer is smaller for disks with larger particle sizes. As the particles grow larger, they sediment more to the midplane, and the surface layer becomes thinner. On the other hand, for the mixed-sized particles (Fig. 4$`b`$), the mass of the surface layer does not vary significantly with the maximum size. Because smaller particles exist at higher altitudes above the midplane, the surface layer is composed mainly of the smallest particles, and the quantity of smallest particles in the surface layer does not vary with the maximum size.

Although Figure 3 implies that the outflow velocity at the surface layer is high, the outflow mass flux is not large because of the small mass contained in the surface layer. The disk evolution time $`t_{\mathrm{out}}`$ through the surface outflow is shown in Figure 5$`a`$. The evolution time of the disk composed of $`0.1\mathrm{\mu m}`$ particles is about $`10^7\mathrm{yr}`$ at 1 AU and $`10^8\mathrm{yr}`$ at 100 AU. Thus, in $`10^7`$ yr, the dust outflow mass which has passed 1 AU becomes comparable to the dust disk mass inside 1 AU. If the radial motion of the dust disk occurs only through surface outflow, the particles in the innermost part of the disk will escape from the region inside 1 AU in $`10^7`$ yr, and a hole will appear. The hole expands, and in $`10^8`$ yr, the radius of the hole reaches 100 AU. There is, however, an accretion flow under the surface layer. Inhibition of disk clearing by the accretion flow is discussed in §3.3 below. As the particle size increases, the evolution time increases as well. This is because particles with larger sizes concentrate more at the midplane and the surface layer becomes thinner. The blow-off time $`t_b`$ for a larger particle at its surface layer is also larger because of the higher local density of the gas at the lower surface layer (Fig. 3). If the disk is composed of $`10\mathrm{\mu m}`$ particles, it takes about $`10^9\mathrm{yr}`$ for a disk of radius 100 AU to evolve.

Figure 5$`b`$ shows the evolution times for disks composed of mixed-sized particles. The disk containing $`[0.01\mathrm{\mu m}1\mathrm{\mu m}]`$ particles evolves in about $`3\times 10^7`$ yr at 1 AU and in $`3\times 10^8`$ yr at 100 AU. The disk evolution occurs mainly through the outflow of $`0.1\mathrm{\mu m}`$ particles. While the surface layer is composed primarily of the smallest ($`0.01\mathrm{\mu m}`$) particles, these particles do not receive radiation pressure as effectively as particles with $`0.1\mathrm{\mu m}`$ (see Fig. 2). Thus, the surface outflow is mainly carried by a minor fraction of particles in the $`0.1\mathrm{\mu m}`$ range, and the disk evolution time of the mixed-sized particles is much longer than that of the disk with the single-sized particles of $`0.1\mathrm{\mu m}`$. The evolution time increases with the maximum particle size. The dependence of the maximum size on the evolution time, however, is not as significant as in the case of single-sized particles. For disks composed of mixed-sized particles, particles of size $`0.1\mathrm{\mu m}`$ form the main contribution to the surface outflow. The number density of $`0.1\mathrm{\mu m}`$ particles decreases as the maximum particle size increases. However, because the thickness of the surface layer is determined by the optical depth from the star and the opacity comes mainly from $`0.1\mathrm{\mu m}`$ particles, the total mass of $`0.1\mathrm{\mu m}`$ particles in the surface layer does not change significantly (Fig. 4$`b`$). Hence, the outflow flux maintains a relatively constant value.

### 3.3 Comparison with the Inward Migration Time

The surface outflow leads to disk clearing only if it dominates the under-surface inflow. The disk evolution time-scale, $`t_{\mathrm{in}}`$, through the inflow is shown as the dashed lines in Figure 5. As discussed in Paper I, for sufficiently small particles that are well coupled to the gas, the inflow velocity is similar to the accretion velocity of the gas. In the standard model, particles smaller than $`1\mathrm{\mu m}`$ have velocities nearly indistinguishable from the gas throughout the disk. Particles of $`10\mathrm{\mu m}`$ decouple from the gas motion at large distances ($`r50`$ AU) and have a shorter evolution time than the gas and smaller particles. The inflowing mass flux is much larger (10 times or more) than the mass flux of the surface outflow even for $`0.1\mathrm{\mu m}`$ particles. Therefore, the net mass flux of particles is dominated by the inflow beneath the surface layer. In the disks composed of mixed-sized particles (Fig. 5$`b`$), the inflow also dominates the net flow. The surface outflow has a very minor effect on the disk evolution in our standard model, regardless of the particle size.

### 3.4 Various Models

In our standard model, the surface outflow is negligible for disk evolution compared to the inflow inside the disk. This flow outcome is due to the surface layer being very thin such that its mass is a small fraction of the entire dust disk mass. The surface outflow may be important under some conditions where the surface layer is thicker or the outflow is much faster than the standard model. In this subsection, we consider several models to see under which conditions the surface outflow is important.

#### 3.4.1 High Luminosity

Massive young stars have a larger luminosity-mass ratio. For example, $`2.5M_{\mathrm{}}`$ stars have about $`20L_{\mathrm{}}/M_{\mathrm{}}`$ for $`10^7`$ yr (D’Antona & Mazzitelli 1994). Under such a high luminosity-mass ratio, the dust particles are exposed to a stronger radiation pressure ($`\beta `$ is proportional to $`L/M`$). Figure 6$`a`$ shows the evolution time of the disk in the case where the stellar luminosity-mass ratio is 15 times larger than the standard model. The outflow velocity is 15 times faster than in the standard model, causing the outflow mass flux to dominate the inflow at large distances ($`r200`$ AU). In this model, dust particles which are initially outside $`200`$ AU are expelled to larger distances, while particles inside $`200`$ AU are accreted onto the star. It is expected that a gap around $`200`$ AU forms in $`10^7`$ yr.

#### 3.4.2 Small Disk Mass

If the total dust mass is smaller, the disk optical depth is smaller so that the layer above the illuminated surface is thicker. Figure 6$`b`$ shows the evolution time of the disks whose masses are $`0.1`$ and $`0.01`$ times smaller than the standard disk. The mass flux of the surface outflow is almost independent of the disk mass because the mass of the surface layer, which is determined by the optical thickness, is approximately constant. Therefore, the relative importance of the surface outflow, which is proportional to the mass ratio of the surface layer to the whole disk, increases with decreasing the total dust mass. The evolution time through the inflow (see the dashed lines) at the inner part of the disk is similar to that of the standard model, while at the outer part it decreases with decreasing disk mass. As the disk mass decreases, the gas density at the outer region of the disk becomes insufficient to keep dust particles moving together with the gas. The inward drift velocity of particles relative to the gas becomes faster for disks with smaller gas mass. For such drifting particles, the evolution time is proportional to $`r^{2+pq/2}`$ (see Paper I, §3.3.3), which is constant with $`r`$ in our models ($`p=2.25`$ and $`q=0.5`$). For the disk with one tenth the mass as the standard model, the evolution time through the surface outflow is still longer than that of the inflow. If the disk is $`10^2`$ times smaller than the standard model (the gas mass inside 100 AU is $`2.5\times 10^4M_{\mathrm{}}=83M_{\mathrm{}}`$), the surface outflow flux is comparable to the inflow flux at around 20 AU. Thus, there is no net flow at radii in the neighborhood of 20 AU. The dust particles initially outside 20 AU move inward until their migration is terminated at 20 AU. At the same time, the particles initially inside 20 AU continue to move inward and accrete onto the star. Therefore, the particle migration induces accumulation of dust at 20 AU and forms a ring there. Gas disks around young stars are expected to dissipate in $`10^7`$ yr (Hartmann et al. 1998). When the gas mass is reduced to about $`100M_{\mathrm{}}`$, a ring appears at a radius of a few tens of AU.

#### 3.4.3 High Diffusivity of Dust

In the standard model, the diffusivity of dust particles subject to turbulent motion of the gas is assumed to be the same as the diffusivity of the gas, i.e., the Schmidt number, Sc, is unity. Experiments show that Sc can be as small as 0.1 and dust particles are much more diffusive than the gas, if the Stokes number, $`\mathrm{St}`$, is of order $`10^2`$ (e.g., Fig. 1 in Cuzzi et al. 1993). The Stokes number is considered to be of the same order as the non-dimensional stopping time $`T_s`$. At radii of several hundred AU, $`T_s10^2`$ for $`0.1\mathrm{\mu m}`$ particles at the illuminated surface ($`z3h_g`$). Thus, it is expected that such particles will spread to higher altitude, where they are easily blown off. Figure 6$`c`$ shows the evolution time for such particles with $`\mathrm{Sc}=0.1`$. The dust particles are more diffusive and are stirred up to higher altitude by turbulence than in the case of $`\mathrm{Sc}=1`$. The illuminated surface is located higher, where the gas density is lower and the outflow velocity is faster. The outflow mass flux, however, is still smaller than the inflow mass flux within $`r500`$ AU. Inflow under the surface layer dominates even for the particles with high diffusivity.

#### 3.4.4 Particle Growth

In the above discussion, we mainly used models which have a power-law size distribution with a maximum size of $`1\mathrm{\mu m}`$. The dust particles in circumstellar disks may have evolved to sizes larger than $`1\mathrm{\mu m}`$ (Beckwith & Sargent 1991; Mannings & Emerson 1994; D’Alessio, Calvet, & Hartmann 2001). Figure 5$`b`$ shows the evolution time of the disks as the maximum size of the particles is varied. The evolution time through the surface outflow does not change significantly as discussed in §3.1. The evolution time through inflow inside the disk, however, rapidly decreases with the growth of particles. At the outer part of the disk where the gas density is insufficient to keep the largest particles moving together with the gas, the inflow velocity of the particles is proportional to their size. For the size distribution with a power-law index $`3.5`$, the dust mass is dominated by the largest particles, and the inflow mass flux is proportional to the size of the largest particle. The evolution time through the inflow is as small as $`10^5`$ yr when the largest particle is $`1`$ mm, which is much shorter than that through the surface outflow. As the particles grow larger, the surface outflow becomes less important for the overall disk evolution.

## 4 Discussion and Summary

### 4.1 Progress in the Disk Evolution

We compared the surface outflow flux to the under-surface inflow flux and discussed conditions for disk clearing. In those comparisons, we assumed the dust density profiles described in §2, and neglected time variations of the dust density. In real disks, however, the dust density varies as the surface outflow proceeds. We discuss how density changes affect the subsequent evolution.

As dust particles above the illuminated surface are blown off by the stellar radiation, new particles continually diffuse into and replenish this region from below. The diffusion time over which the surface layer is re-populate is

$$t_{\mathrm{rp}}\frac{\mathrm{Sc}(3h_g)^2}{\nu }\frac{9}{\alpha }\mathrm{\Omega }^110^3r_{\mathrm{AU}}^{3/2}\mathrm{yr},$$
(35)
where we use the fact that the distance between the surface layer and the midplane is about three disk scale heights for $`\alpha =10^3`$ and $`\mathrm{Sc}=1`$. The diffusion time ($`10^6`$ yr at 100 AU) is much shorter than the disk evolution time $`t_{\mathrm{out}}`$. Thus during the disk evolution, the dust in the surface layer is continuously replenished, and the surface outflow continues until the dust in the midplane is depleted.

In disks composed of mixed-sized particles, small ($`s0.1\mathrm{\mu m}`$) particles constitute the majorities of the surface outflow. The small particles are rapidly removed from the disk, and are continuously replenished through the collisional destruction of larger particles. The collisional time scale, $`t_{\mathrm{col}}=(\pi s^2n_dv_{\mathrm{rel}})^1`$, where $`n_d`$ is the particle number density and $`v_{\mathrm{rel}}`$ is the relative velocity of particles, is much shorter than the disk evolution time. For example, if we assume that the majority of the dust is composed of 1 cm particles, the relative velocity is about $`10\mathrm{cm}\mathrm{s}^1`$ at 1 AU (from Fig. 3 in Weidenschilling & Cuzzi 1993), and the collisional time is $`t_{\mathrm{col}}10^3`$ yr. If the dust is composed of smaller particles, the collisional time is smaller. The outcome of the collisions, however, is not well understood, and we do not know how many small particles are produced through the collisional destruction. In this paper, we simply assume that the particle size distribution is in an equilibrium state. The minimum size, the maximum size, and the power-law index of the distribution do not vary with time. In a realistic situation, however, at least the maximum size would grow with time, leading to inhibition of the disk clearing, as discussed in §3.4.4. Calculations that take account of dust growth are needed to solve this issue.

Because the disk evolution time $`t_{\mathrm{out}}`$ is shorter at smaller radii, the dust outflow proceeds from the inner disk. The particles blown from the inner disk accumulate in the outer disk. This process increases the dust-to-gas ratio at large radii. We can ask whether this increase in the dust-to-gas ratio affects the subsequent disk evolution. Suppose that the dust particles inside a radius $`r`$ are removed and an inner hole (or a gap) is produced. The removed particles pile up in an annulus just outside the inner hole. The width of the annulus is larger than the disk scale height. On the diffusion time $`t_{\mathrm{rp}}`$, the accumulated particles diffuse over a few disk scale heights both in the vertical and in the radial directions, and during a disk evolution time $`t_{\mathrm{out}}>t_{\mathrm{rp}}`$, the particles diffuse even further. Thus, the surface density of the annulus increases at most by a factor $`\pi r^2\mathrm{\Sigma }_{d,\mathrm{all}}/(2\pi r\mathrm{\Sigma }_{d,\mathrm{all}}dr)r/h_g10`$, where $`drh_g0.1r`$ is the width of the annulus. Because the gas density does not change, the dust-to-gas ratio of the annulus becomes $`10`$ times larger. The opacity $`\kappa `$, which is defined for a unit mass of the gas and is proportional to the dust-to-gas ratio, also increases by a factor of order $`10`$ times. (If the particles diffuse as far as the inner hole radius $`r`$ in a disk evolution time $`t_{\mathrm{out}}`$, then the hole is refilled and dust removal is inhibited. The time-scale for the refill is comparable to the gas accretion time-scale $`r^2/\nu `$, which is shown as the dotted lines in Figure 5. In the standard model, the refill effectively inhibits the dust removal. Only in the high stellar luminosity or the low mass disk cases, the surface outflow can dominate the refill.)

Figure 7 shows how the increase in the dust-to-gas ratio affects the outflow mass flux. In Figure 7, the outflow flux in the cases for $`f_{\mathrm{dust}}=0.1`$ (i.e., the dust-to-gas ratio is 10 times larger than the standard model) and for $`f_{\mathrm{dust}}=10^3`$ (10 times smaller dust-to-gas ratio) are shown with the standard model. It is seen that the outflow flux depends only weakly on the dust-to-gas ratio. The change in the dust-to-gas ratio has a slight effect on the location of the surface layer, but it hardly changes the dust mass of the surface layer. The thickness of the surface layer is determined so as to make the optical depth be unity, and its mass content does not depend on the dust-to-gas ratio. The variation in the outflow flux arises mainly from change in the location of the surface layer and from change in the angle of the incident starlight to the surface layer, and these effects are small. Therefore, the accumulation of blown off particles increases the dust-to-gas ratio at the disk edge just outside the inner hole, but does not significantly decrease the outflow flux there. Even after the increase in the dust-to-gas ratio in the outer disk, the surface outflow proceeds and the disk evolution continues with the time-scale plotted in Figures 5 and 6.

### 4.2 Gap Formation in Herbig Ae/Be Disks

The phenomenon of surface stripping can effectively evolve the disk structure around young luminous stars whose luminosity-mass ratios are larger than $`15`$ times the solar value, as shown in §3.4.1. In the outer part ($`200`$ AU) of such disks, the surface outflow dominates the inflow through the midplane, causing a net outward flow. On the other hand, within 200 AU, the dust flows inward onto the central star. Around 200 AU, the dust density decreases. The outflow or inflow flux around 200 AU is strong enough that it can carry the total amount of the dust inside 200 AU in $`10^7`$ yr. This flux maintains a relatively constant value throughout the evolution even after the dust density and the dust-to-gas ratio have decreased (see Fig. 7). Consequently, the dust density declines at an almost constant rate. The time-scale for dust diffusion to transport dust from the midplane to replenish the surface layer is about $`3\times 10^6`$ yr at 200 AU. We expect the dust density at 200 AU will decline considerably in $`10^7`$ yr.

The formation of gaps can be examined by radio observations. At 200 AU, the gas surface density in the standard model is $`\mathrm{\Sigma }_g=1.8\mathrm{g}\mathrm{cm}^2`$, which is optically thin in the vertical direction in radio continuum region around 1 mm ($`\tau _{1\mathrm{m}\mathrm{m}}10^2`$). Therefore, the density decline will be observed as a decrease in the optical depth.

The disk is highly optically thick in the visible or near infrared light ($`\tau _{\mathrm{vis}}200`$ at 200 AU). Thus, the density decrease is difficult to examine directly using visible or near-infrared observations. If the gap region is shadowed by the inner disk, or if the edge of the outer disk is illuminated by the central star, the variation in the surface irradiation will probably be observed. There are a few Herbig Ae/Be stars that show gaps in the scattered light images of their disks. The disk around HD 163296 has a dark lane at 325 AU (Grady et al. 2000). HD 100546 also has a disk with a gap at 250 AU (Grady et al. 2001). In order to make direct comparisons with these observations, we need further modeling that takes account of the time evolution of the dust density and that incorporates simulated scattered light images.

In this paper, we discussed gap formation through the surface stripping caused by radiation pressure. It is useful to discuss how to distinguish this mechanism from the other gap formation processes such as the planetary perturbations (Lin & Papaloizou 1986; Takeuchi, Miyama & Lin 1996). Obviously, the surface stripping itself cannot generate non-axisymmetric features, whereas an embedded planet can. In gas-free disks, a planet forms characteristic disk structures that arises from orbital resonances (Roques et al. 1994; Dermott et al. 1994; Liou & Zook 1999; Wyatt et al. 1999; Ozernoy et al. 2000; Wilner et al. 2002; Quillen & Thorndike 2002), and in gas disks the planet excites spiral density waves (Miki 1982; Sekiya, Miyama, & Hayashi 1987; Bryden et al. 1999; Kley 1999; Lubow, Seibert, & Artymowicz 1999; Miyoshi et al. 1999; Tanigawa & Watanabe 2002). The detection of such non-axisymmetric features would rule out surface stripping as a gap formation mechanism. The other notable feature of surface stripping is that it does not affect the gas density structure. The gas disks still have monotonic density structure even after the dust disks have formed gaps. Therefore, detection of gaps in gas disks also rules out surface stripping. In order to explain the cause of the gaps around HD 163296 and HD 100546, it is important to determine with future observations whether their gas disks have gaps or not.

### 4.3 Summary

Dust particles that are directly irradiated by the star move outward under the combined action of gas drag and radiation pressure. In this paper, the surface outflow of particles in optically thick disks is studied. Our results are as follows.

1. Particles in the surface layer, which are directly exposed to the starlight, move outward. In a disk resembling the minimum-mass-solar-nebula, the height of the surface layer is about three scale heights of the gas disk and the mass is $`10^6`$ times the total dust mass at 1 AU and $`10^4`$ times the total dust mass at 100 AU. For this case, the mass flux of the surface outflow is much smaller than the inflow mass flux inside the disk, and the surface outflow does not contribute significantly to the dust disk evolution.

2. The surface outflow dominates the inflow in the following cases: ($`a`$) If the luminosity-mass ratio of the central star is 15 times larger than the solar value, the velocity of the surface outflow is large enough at the outer part of the disk ($`r200`$ AU). ($`b`$) If the mass of the gas disk is as small as $`10^2`$ times the minimum-mass disk ($`100M_{\mathrm{}}`$ inside 100 AU), the mass ratio of the surface layer to the whole dust disk is large enough at a few tens AU. In such cases, the surface outflow induces formation of gaps or rings.

3. The surface outflow is mainly composed of particles with size of order $`0.1\mathrm{\mu m}`$. As the particles grow, the largest particles around the midplane move inward more rapidly, while the outflow flux at the surface remains at a near-constant value. This results in a diminishing importance of the surface outflow for the overall disk evolution.

We would like to thank the referee for helpful comments. We also thank Greg Laughlin for his careful reading the manuscript and correcting grammar. This work was supported in part by an NSF grant AST 99 87417 and in part by a special NASA astrophysical theory program that supports a joint Center for Star Formation Studies at UC Berkeley, NASA-Ames Research Center, and UC Santa Cruz. This work is also supported by NASA NAG5-10612 through its Origin program and by JPL 1228184 through its SIM program.


format: text

text:
# Adiabatic quantum search algorithm for structured problems

## I Introduction

Grover’s quantum algorithm solves an unstructured search problem in a time of order $`\sqrt{N}`$, where $`N`$ is the dimension of the search space, which corresponds to a quadratic speed-up over a classical search Grover (1996). This algorithm is proved to be optimal in the case of unstructured search problems Zalka (1999). Naturally, it can also be used to solve a structured search problem with a quadratic speed-up over a naive classical search that would exhaustively check every possible solution. However, exploiting the structure of the problem is well known to lead to better classical search algorithms. It is therefore tempting to imagine that better quantum search algorithms may be devised similarly by exploiting the problem structure. Following this, Cerf, Grover and Williams have shown that this could be done by partitioning the unknown variables into two (or more) sets and nesting a quantum search over one set into another search over two (or more) sets, yielding an average complexity of order $`\sqrt{N^\alpha }`$, with $`\alpha <1`$ J.Cerf et al. (2000).

While this algorithm, as well as Grover’s original algorithm, stay within the standard paradigm of quantum computation based on quantum circuits, a new type of quantum algorithms based on adiabatic evolution has been introduced lately Farhi et al. (2000); van Dam et al. (2001). In particular, a quantum adiabatic analogue of Grover’s search algorithm has been proposed in Ref. Roland and Cerf (2002), which works for unstructured search problems. The use of quantum adiabatic algorithms has also been analyzed for solving structured problems such as $`k`$-SAT, but in such a way that until now only a numerical study has been possible Farhi et al. (2001). Recently, the study of quantum adiabatic algorithms has progressed even further after Aharonov and Ta-Shma demonstrated their universality, namely that any problem that may be solved efficiently in the standard (circuit-based) model of quantum computation can also be solved efficiently by a quantum adiabatic algorithm Aharonov and Ta-Shma (2003). This provides a strong incentive to the search for new quantum adiabatic algorithms.

The purpose of this paper is to bring the ideas of nested quantum search and quantum adiabatic computation together, in order to devise a new quantum adiabatic algorithm adapted to structured problems. More specifically, we will show that an adiabatic search over a subset of the variables can be used to build a better initial Hamiltonian for the global adiabatic search. With this adiabatic algorithm, we recover the same complexity as the nested circuit-based algorithm of Ref. J.Cerf et al. (2000).

## II Adiabatic theorem

Let us briefly recall the Adiabatic Theorem and how it may be used to design quantum algorithms by adiabatic evolution.

We know that if a quantum system is prepared in the ground state of a time-independent driving Hamiltonian, it remains in this state. The Adiabatic Theorem states that if this Hamiltonian becomes time-dependent, the system will still stay close to its instantaneous ground state as long as the variation is slow enough.

More specifically, if $`|E_0;t`$ and $`|E_1;t`$ are the ground and first excited states of the Hamiltonian $`H(t)`$, with energies $`E_0(t)`$ and $`E_1(t)`$, we define the minimum gap between these eigenvalues

$$g_{\mathrm{min}}=\underset{0tT}{\mathrm{min}}\left[E_1(t)E_0(t)\right]$$
(1)
and the maximum value of the matrix element of $`dH/dt`$ between the eigenstates as

$$D_{\mathrm{max}}=\underset{0tT}{\mathrm{max}}\left|\frac{dH}{dt}_{1,0}\right|$$
(2)
with $`\frac{dH}{dt}_{1,0}=E_1;t|\frac{dH}{dt}|E_0;t`$. The Adiabatic Theorem states that if we prepare the system at time $`t=0`$ in its ground state $`|E_0;0`$ and let it evolve under the Hamiltonian $`H(t)`$, then

$$|E_0;T|\psi (T)|^21\epsilon ^2$$
(3)
provided that

$$\frac{D_{\mathrm{max}}}{g_{\mathrm{min}}^2}\epsilon ,$$
(4)
where $`\epsilon 1`$.

Now, we may apply to the system a Hamiltonian for which the ground state encodes the unknown solution of a problem. According to the Adiabatic Theorem, we know that we may get the system very close to this solution state by preparing it in the (known) ground state of another Hamiltonian, and then by progressively changing it to the Hamiltonian of our problem. This simple idea is central to the quantum algorithms by adiabatic evolution Farhi et al. (2000); van Dam et al. (2001).

## III Quantum search by local adiabatic evolution

As exposed in Roland and Cerf (2002), this principle may be used to perform a quantum search. Suppose that among $`N`$ states, we have to find the $`M`$-times degenerate ground state of a Hamiltonian

$$H_f=I\underset{m}{}|mm|,$$
(5)
where $``$ is the ensemble of solutions (of size $`M`$). We initially prepare the system in an equal superposition of all could-be solutions:

$$|s=\frac{1}{\sqrt{N}}\underset{i𝒩}{}|i.$$
(6)
This superposition is the ground state of the following Hamiltonian:

$$H_i=I|ss|.$$
(7)
We now apply $`H_i`$ to the system and switch adiabatically to $`H_f`$. If we perform an adiabatic evolution

$$H(t)=(1s(t))H_i+s(t)H_f,$$
(8)
where $`s(t)`$ is a (carefully chosen) monotonic function with $`s(0)=0`$ and $`s(T)=1`$, we will finally obtain a state close to a ground state of $`H_f`$:

$$|\psi _f\frac{1}{\sqrt{M}}\underset{m}{}|m$$
(9)
as long as

$$T=O\left(\sqrt{\frac{N}{M}}\right)$$
(10)
This algorithm is refered to as local because $`s(t)`$ is chosen such that the Adiabatic Theorem is obeyed locally, at each time (see Roland and Cerf (2002) for details).

## IV Structured problems

In this article, we consider a class of problems where one has to find an assignment for a set of variables. For each additional variable considered, new constraints appear and reduce the set of satisfying assignments. This corresponds to most problems encountered in practice ($`k`$-SAT, graph coloring, planning, combinatorial optimization, …).

For a set of $`n_A`$ variables denoted as $`A`$, there is a corresponding set of constraints $`C_A`$. We may define a function $`f_A`$ that tells if an assignment of the variables in $`A`$ satisfies the constraints in $`C_A`$.

$`f_A`$ $`:`$ $`(_d)^{n_A}\{0,1\}`$ (13)
$`:`$ $`x\{\begin{array}{cc}0\hfill & \text{if }x\text{ does not satisfy }C_A\hfill \\ 1\hfill & \text{if }x\text{ satisfies }C_A\text{,}\hfill \end{array}`$
where $`d`$ is the number of possible assignments for each variable ($`d=2`$ for bits). As quantum gates have to be reversible, the quantum equivalent of this function will be an oracle:

$$O_A:_{N_A}_2_{N_A}_2:|x|y|x|yf_A(x),$$
(14)
where $`N_A=d^{n_A}`$. It is shown in Ref. Roland and Cerf (2003) that this oracle is closely related to a Hamiltonian whose ground states, of energy $`0`$, are the basis states encoding a satisfying assignment and whose excited states, of energy $`1`$, are all other basis states:

$$H_A|x=\{\begin{array}{cc}|x\hfill & \text{if }x\text{ does not satisfy }C_A\hfill \\ 0\hfill & \text{if }x\text{ satisfies }C_A\hfill \end{array}$$
(15)
or

$$H_A=I_A\underset{m_A_A}{}|m_Am_A|,$$
(16)
where $`_A`$ is the set of satisfying assignments for the variables in $`A`$. It is possible to reproduce the application of this Hamiltonian during a certain time by using two calls to the quantum oracle $`O_A`$ (see Roland and Cerf (2003) for details).

Now suppose we consider a larger set of variables $`n_{AB}=n_A+n_B`$ that have to satisfy a set of constraints $`C_{AB}C_A`$. To discriminate between assignments satisfying $`C_{AB}`$ or not, we will use an oracle $`O_{AB}`$ or a corresponding Hamiltonian $`H_{AB}`$ defined as in Eqs. (14) and (15). The basic idea of our structured search will be to find the solutions of $`C_{AB}`$ by first building the assignments of the $`n_A`$ primary variables satisfying $`C_A`$, then by completing them with all possible assignments of the $`n_B`$ secondary variables, and finally by searching among these could-be solutions the global satisfying assignments.

## V Structured search with two levels of nesting

This problem is of the same type as the one considered in J.Cerf et al. (2000), for which the technique of nesting was introduced in the context of the traditional implementation of Grover’s algorithm on a conventional quantum circuit. Here, we apply this technique to the adiabatic quantum search algorithm.

Suppose we divide the variables of our problem into two subsets $`A`$ ($`n_A`$ elements) and $`B`$ ($`n_B`$ elements). First, we will perform a search on the variables in $`A`$ using the Hamiltonian $`H_A`$ that encodes the constraints in $`C_A`$:

$$H_A=I_A\underset{m_A_A}{}|m_Am_A|.$$
(17)
Then we will use the Hamiltonian $`H_{AB}`$ acting on all variables in $`AB`$ and encoding the whole set of constraints $`C_{AB}`$

$`H_{AB}`$ $`=`$ $`I_{AB}`$
$`{\displaystyle \underset{(m_A,m_B)_{AB}}{}}|m_Am_A||m_Bm_B|`$
to construct a superposition of the solutions of the full problem $`_{AB}`$. A final measurement of the quantum register then gives one of the global solutions at random.

### V.1 Adiabatic search on the primary variables

The preliminary search on the variables in $`A`$ is a simple unstructured search as explained in section III. As there are $`n_A`$ variables in $`A`$, the corresponding Hilbert space is of dimension $`N_A=d^{n_A}`$. Let $`M_A`$ be the number of solutions in $`_A`$. Performing an adiabatic quantum search, we may thus transform the initial state

$$|s_A=\frac{1}{\sqrt{N_A}}\underset{i𝒩_A}{}|i_A$$
(19)
into a state close to the uniform superposition of all solutions in $`_A`$

$$|\psi _{m_A}=\frac{1}{\sqrt{M_A}}\underset{m_A_A}{}|m_A$$
(20)
in a time of order

$$T_A=0\left(\sqrt{\frac{N_A}{M_A}}\right).$$
(21)
### V.2 Adiabatic search on the secondary variables

We will now perform a preliminary search in the Hilbert space of dimension $`N_B=d^{n_B}`$ of the secondary variables in $`B`$ by extending the partial solutions $`|m_A`$. We prepare the variables in $`B`$ in a state that is the uniform superposition

$$|s_B=\frac{1}{\sqrt{N_B}}\underset{j𝒩_B}{}|j_B.$$
(22)
Globally, the system is thus in the superposition:

$`|\psi _0_{AB}`$ $`=`$ $`|\psi _{m_A}|s_B`$ (23)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{M_AN_B}}}{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{m_A_A}{j𝒩_B}}{}}|m_A|j_B,`$
where some terms correspond to a global solution of the problem \[$`(m_A,j)_{AB}`$ satisfying all constraints in $`C_{AB}`$\] and the others to a partial solution only \[$`m_A_A`$ satisfies $`C_A`$ but $`(m_A,j)_{AB}`$ does not satisfy $`C_{AB}`$\]. We now divide the set $`_A`$ of solutions of $`C_A`$ into two subsets: $`_A^S`$ will be the set of $`m_A`$’s for which there exists at least one solution $`(m_A,m_B)`$ of $`C_{AB}`$ and $`_A^{NS}`$ the set of $`m_A`$’s for which there is no such solution.

$`_A^S`$ $`=`$ $`\{m_A_A|m_B,(m_A,m_B)_{AB}\}`$ (24)
$`_A^{NS}`$ $`=`$ $`\{m_A_A|j,(m_A,j)_{AB}\}`$ (25)
Of course, we thus have $`_A=_A^S_A^{NS}`$. We may now rewrite our initial state (23) as

$`|\psi _0_{AB}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{M_AN_B}}}{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{m_A_A^{NS}}{j𝒩_B}}{}}|m_A|j_B`$ (26)
$`+{\displaystyle \frac{1}{\sqrt{M_AN_B}}}{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{m_A_A^S}{j𝒩_B}}{}}|m_A|j_B.`$
In the first part of this expression, no term correspond to a solution of the full problem, whereas in the second part, some terms do and others do not. The goal of this stage of the computation will be to increase the amplitude of the solution terms in this last part. To achieve this, we perform an adiabatic evolution using as initial Hamiltonian

$$H_i=I_A(I_B|s_Bs_B|),$$
(27)
that has $`|\psi _0_{AB}`$ as a ground state. The final Hamiltonian will be

$$H_f=H_{AB}H_AI_B.$$
(28)
We see that these Hamiltonians share the following properties:

1. They do not act on states $`|i_A|s_B`$ corresponding to assignments $`i`$ of $`𝒩_A`$ that do not satisfy $`C_A`$:
$`H_{i,f}|i_A|s_B=0i_A`$.
2. They do not couple states corresponding to different $`m_A`$’s:
$`{}_{B}{}^{}j|m_A|H_{i,f}|m_A^{}|j^{}_{B}^{}=0m_Am_A^{}_A,j,j^{}𝒩_B`$.

It follows that the instantaneous Hamiltonian of the adiabatic evolution $`H(t)`$ satisfies the same properties. Keeping this in mind, it may easily be shown that the effect of the adiabatic evolution will be to perform independent adiabatic searches for each $`m_A_A`$. More precisely, each term in $`|\psi _0_{AB}`$

$$\frac{1}{\sqrt{N_B}}\underset{j𝒩_B}{}|m_A|j_B$$
(29)
will evolve to

$$\frac{1}{\sqrt{M_{B/m_A}}}\underset{m_B_{B/m_A}}{}|m_A|m_B,$$
(30)
as long as

$$T_{m_A}=O\left(\sqrt{\frac{N_B}{M_{B/m_A}}}\right),$$
(31)
where $`_{B/m_A}`$ is the set of $`m_B`$’s such that $`(m_A,m_B)_{AB}`$ and $`M_{B/m_A}`$ is the number of these elements. For this condition to be satisfied for all $`m_A`$’s simultaneously, we must take

$$T_B=\underset{m_A}{\mathrm{max}}T_{m_A}=O\left(\sqrt{\frac{N_B}{\mathrm{min}_{m_A}M_{B/m_A}}}\right).$$
(32)
At the end of this second stage, we thus have constructed a state close to

$`|\psi _{AB}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{M_AN_B}}}{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{m_A_A^{NS}}{j𝒩_B}}{}}|m_A|j_B`$ (33)
$`+{\displaystyle \frac{1}{\sqrt{M_A}}}{\displaystyle \underset{m_A_A^S}{}}e^{i\varphi _{m_A}}|m_A`$
$`{\displaystyle \frac{1}{\sqrt{M_{B/m_A}}}}{\displaystyle \underset{m_B_{B/m_A}}{}}|m_B`$
$`=`$ $`\sqrt{{\displaystyle \frac{M_A^{NS}}{M_A}}}|\psi ^{NS}+\sqrt{{\displaystyle \frac{M_A^S}{M_A}}}|\psi ^S,`$
where $`\varphi _{m_A}`$’s are phases appearing during the evolution,

$`|\psi ^{NS}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{M_A^{NS}N_B}}}{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{m_A_A^{NS}}{j𝒩_B}}{}}|m_A|j_B`$ (34)
$`|\psi ^S`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{M_A^S}}}{\displaystyle \underset{m_A_A^S}{}}e^{i\varphi _{m_A}}|m_A`$ (35)
$`{\displaystyle \frac{1}{\sqrt{M_{B/m_A}}}}{\displaystyle \underset{m_B_{B/m_A}}{}}|m_B`$
and $`M_A^{NS}`$ (resp. $`M_A^S`$) is the number of elements in set $`_A^{NS}`$ (resp. $`_A^S`$).

### V.3 Global adiabatic search

Stages A and B define a unitary evolution $`U`$ that applies the initial state $`|s_A|s_B`$ onto $`|\psi _{AB}`$:

$`U|s_A|s_B`$ $``$ $`|\psi _{AB}`$ (36)
$`=`$ $`\sqrt{{\displaystyle \frac{M_A^{NS}}{M_A}}}|\psi ^{NS}+\sqrt{{\displaystyle \frac{M_A^S}{M_A}}}|\psi ^S.`$ (37)
In this state, we now need to decrease the amplitude of the first term, corresponding to partial solutions only, and increase the amplitude of the second term, corresponding to global solutions. This could be realized efficiently by performing an adiabatic search using as initial Hamiltonian:

$`H_i`$ $`=`$ $`I_{AB}|\psi _{AB}\psi _{AB}|`$ (38)
$``$ $`U(I_{AB}|s_As_A||s_Bs_B|)U^{}`$ (39)
$``$ $`UH_0U^{},`$ (40)
where $`H_0=I_{AB}|s_As_A||s_Bs_B|`$, and as final Hamiltonian

$`H_f`$ $`=`$ $`H_{AB}`$
$`=`$ $`I_{AB}{\displaystyle \underset{(m_A,m_B)_{AB}}{}}|m_Am_A||m_Bm_B|`$
during a time

$$T_C=0\left(\sqrt{\frac{M_A}{M_A^S}}\right).$$
(42)
Unfortunately, we do not have access to $`H_i`$, so that the interpolating Hamiltonian $`H(s)=(1s)H_i+sH_f`$ cannot be applied directly. However, we will see in Section VI that the basic steps of the quantum circuit implementation of this adiabatic algorithm only require the application of $`H_i`$ during a particular time $`t`$, that is

$$e^{iH_it}e^{iUH_0U^{}t}=Ue^{iH_0t}U^{}.$$
(43)
Hence, each application of $`H_i`$ during a time $`t`$ will be equivalent to sequentially applying $`U^{}`$, $`e^{iH_0t}`$, and $`U`$, which means performing the adiabatic evolution $`U`$ (stages A and B) backwards, then applying $`H_0`$ for a time $`t`$, and finally rerun $`U`$ forwards (stages A and B).

In section VI, we will see that, when discretizing the evolution, we must take a number of steps $`r_C`$ of order $`T_C`$. We may now evaluate the complexity of this algorithm. As it consists of $`r_C`$ steps, each involving two applications of $`U`$ or $`U^{}`$, that last a time of order $`T_A+T_B`$, the algorithm finally takes a time of order:

$`T`$ $`=`$ $`(T_A+T_B)r_C`$
$`=`$ $`O\left(\left(\sqrt{{\displaystyle \frac{N_A}{M_A}}}+\sqrt{{\displaystyle \frac{N_B}{\mathrm{min}_{m_A}M_{B/m_A}}}}\right)\sqrt{{\displaystyle \frac{M_A}{M_A^S}}}\right)`$
$`=`$ $`O\left(\sqrt{{\displaystyle \frac{N_A}{M_A^S}}}+\sqrt{{\displaystyle \frac{M_AN_B}{M_A^S\mathrm{min}_{m_A}M_{B/m_A}}}}\right).`$ (45)
Let us notice that with the same hypothesis as in Ref. J.Cerf et al. (2000), namely

$$M_{B/m_A}=1m_A,$$
(46)
then $`M_A^S=M_{AB}`$, and the computation time is

$$T=O\left(\frac{\sqrt{N_A}+\sqrt{M_AN_B}}{\sqrt{M_{AB}}}\right)$$
(47)
i.e., the same complexity as the equivalent circuit-based algorithm exposed in Ref. J.Cerf et al. (2000). A more detailed analysis of this complexity will be performed in section VII.

## VI Discretizing the adiabatic evolution

### VI.1 General method

The method to implement a global adiabatic evolution algorithm on a discrete quantum circuit was initially exposed in van Dam et al. (2001), and was extended to the case of a local adiabatic evolution algorithm in Roland and Cerf (2003). Let us quickly review this method, that uses two successive approximations.

The first approximation consists in cutting the evolution time $`T`$ in $`r`$ intervals $`\mathrm{\Delta }T=\frac{T}{r}`$ and replacing the continuously varying Hamiltonian $`H(t)`$ by a Hamiltonian $`H^{}(t)`$ that is constant during each interval $`\mathrm{\Delta }T`$ and varies at times $`t_j=j\mathrm{\Delta }T`$ only.

$$H^{}(t)=H(t_j)\text{if}t_{j1}tt_j.$$
(48)
It is shown in Roland and Cerf (2003) that for $`H(s)=(1s)H_i+sH_f`$ with $`s=s(t)`$, this approximation introduces a global error for the evolution

$$|U(T)U^{}(T)|_2\sqrt{2\frac{T}{r}|H_iH_f|_2},$$
(49)
where $`|A|_2=\mathrm{max}_{|x=1}A|x`$ is the operator norm of $`A`$. Our algorithm now requires $`r`$ steps of the form

$$U_j^{}=e^{iH(t_j)\mathrm{\Delta }T}=e^{i(1s_j)H_i\mathrm{\Delta }Tis_jH_f\mathrm{\Delta }T},$$
(50)
where $`s_j=s(t_j)`$. As we are able to apply $`H_i`$ and $`H_f`$ separately but not necessarily a simultaneous combination of them, we will approximate $`U_j^{}`$ by

$$U_j^{\prime \prime }=e^{i(1s_j)H_i\mathrm{\Delta }T}e^{is_jH_f\mathrm{\Delta }T}.$$
(51)
This will result in an error

$$|U^{}(T)\underset{j}{}U_j^{\prime \prime }|_2O(\frac{T^2}{r}|[H_i,H_f]|_2)$$
(52)
(see Roland and Cerf (2003) for details).

### VI.2 Application to a structured quantum search

We now consider the case of a structured quantum search. We could apply the discretization procedure to all three stages (A, B, C) of our algorithm in order to implement it on a quantum circuit, but we will concentrate on stage C which is the only one that requires discretization. Nonetheless, it is easy to show that stage A (resp. B) would require a number of steps $`r_A`$ (resp. $`r_B`$) of the same order as the computation time $`T_A`$ (resp. $`T_B`$).

For the final stage, the global adiabatic search, the Hamiltonians $`H_i`$ and $`H_f`$ are defined in Eq. (38)-(V.3). Evaluating the errors introduced by the approximations, we find:

$`|H_iH_f|_2`$ $`<`$ $`1`$ (53)
$`|[H_i,H_f]|_2`$ $`<`$ $`\sqrt{{\displaystyle \frac{M_A^S}{M_A}}}`$ (54)
and, as $`T_C=0\left(\sqrt{\frac{M_A}{M_A^S}}\right)`$,

$`|U(T)U^{}(T)|_2`$ $``$ $`O\left(\sqrt{2{\displaystyle \frac{\sqrt{\frac{M_A}{M_A^S}}}{r}}}\right)`$ (55)
$`|U^{}(T){\displaystyle \underset{j}{}}U_j^{\prime \prime }|_2`$ $``$ $`O\left({\displaystyle \frac{\sqrt{\frac{M_A}{M_A^S}}}{r}}\right).`$ (56)
Therefore, as announced in section V, we have to cut our evolution in a number of steps $`r_C=O\left(\sqrt{\frac{M_A}{M_A^S}}\right)`$ of the same order as $`T_C`$. Each step $`j`$ will take the form:

$`U_j^{\prime \prime }`$ $`=`$ $`e^{i(1s_j)H_i\mathrm{\Delta }T}e^{is_jH_f\mathrm{\Delta }T}`$ (57)
$``$ $`Ue^{i(1s_j)H_0\mathrm{\Delta }T}U^{}e^{is_jH_f\mathrm{\Delta }T},`$ (58)
where the applications of Hamiltonians $`H_0`$ during a time $`(1s_j)\mathrm{\Delta }T`$ and $`H_f`$ during a time $`s_j\mathrm{\Delta }T`$ may be realized by the procedure exposed in Roland and Cerf (2003).

## VII Complexity analysis

To estimate the efficiency of this algorithm, we will follow the same development as in J.Cerf et al. (2000): as we have seen in Section V, under the assumption (46) that we will consider here, the complexity of this adiabatic algorithm has exactly the same form as its circuit-based counterpart.

First of all let us define a few concepts (for details here and throughout this section, we refer the reader to Ref. J.Cerf et al. (2000)). The structured search problem is to find an assignment of $`n_{AB}=n_A+n_B`$ variables among $`d`$ possibilities and satisfying $`e`$ constraints, each involving at most $`k`$ of these variables. We define as a ground instance an assignment of all the variables involved in a particular constraint. A ground instance will be said to be nogood if it violates the constraint. Let $`\xi `$ be the number of those no-good ground instances.

Empirical studies show that the difficulty of solving a structured problem essentially depends on four parameters: the number of variables $`n_{AB}`$, the number of possible assignment per variable $`d`$, the number of variables per constraint $`k`$, and the number of no-good ground instances $`\xi `$. Intuitively, we understand that if $`\xi `$ is small, there are many assignments satisfying the constraints so the problem is easy to solve. On the contrary, if $`\xi `$ is large, the problem is over-constrained and it is easy to show that there is no solution. More precisely, it may be shown that for fixed $`n_{AB}`$ and $`d`$, the average difficulty may be evaluated by the parameter $`\beta =\xi /n_{AB}`$. The hard problems will be concentrated around a critical value $`\beta _c`$.

Let us now estimate the complexity (47). Let $`p(n)`$ be the probability that a randomly generated assignment of the $`n`$ first variables satisfies all the constraints involving these variables. We then have $`M_A=p(n_A)d^{n_A}`$ and $`M_{AB}=p(n_{AB})d^{n_{AB}}`$ while it is shown in J.Cerf et al. (2000) that:

$$p(n)d^{n_{AB}(\beta /\beta _c)(n/n_{AB})^k}.$$
(59)
Eq. (47) becomes

$$T=O\left(\frac{\sqrt{d^{n_A}}+\sqrt{d^{n_{AB}[1(\beta /\beta _c)(n_A/n_{AB})^k]}}}{\sqrt{d^{n_{AB}}(1\beta /\beta _c)}}\right)$$
(60)
or, with $`a=\sqrt{d^{n_{AB}}}`$ and $`x=n_A/n_{AB}`$:

$$T=O\left(\frac{a^x+a^{1\frac{\beta }{\beta _c}x^k}}{a^{1\frac{\beta }{\beta _c}}}\right).$$
(61)
We now optimize $`x`$, the fraction of variables for which we perform a partial search, to minimize the computation time. We have to solve the equation

$$\frac{\beta }{\beta _c}kx^{k1}=a^{\frac{\beta }{\beta _c}x^k+x1}$$
(62)
which, for large $`a`$ (that is large $`n_{AB}`$) approximately reduces to

$$\frac{\beta }{\beta _c}x^k+x1=0.$$
(63)
The solution of this equation $`\alpha `$ ($`0\alpha 1`$) corresponds to the optimal partial search we may perform such that the complexity grows with the smallest power in $`d`$ for $`n_{AB}\mathrm{}`$. This optimal computation time may then be written as

$$T=O\left(\frac{2a^\alpha }{a^{1\frac{\beta }{\beta _c}}}\right)=O\left(\frac{\sqrt{d^{\alpha n_{AB}}}}{\sqrt{d^{n_{AB}(1\frac{\beta }{\beta _c})}}}\right).$$
(64)
Let us now consider the hardest problems for which $`\beta \beta _c`$. For these problems, the complexity reads:

$$T=O\left(\sqrt{d^{\alpha n_{AB}}}\right)$$
(65)
that we may immediately compare to the complexity of an unstructured quantum search $`O\left(\sqrt{d^{n_{AB}}}\right)`$. The gain in the exponent $`\alpha `$ depends on $`k`$ through Eq. (63) For instance, we find $`\alpha =0.62`$ for $`k=2`$, $`\alpha =0.68`$ for $`k=3`$ and $`\alpha 1`$ when $`k\mathrm{}`$.

To compare these results with a classical algorithm, let us consider a specific problem, the satisfiability of boolean formulas in conjunctive normal form, or $`k`$-SAT. For $`3`$-SAT, which is known to be NP-complete, the best known classical algorithm has a worst-case running time that scales as $`O\left(2^{0.45n_{AB}}\right)`$ Paturi et al. (1998), while, as $`\alpha =0.68`$ for $`k=3`$, our quantum adiabatic algorithm has a computation time of order $`O\left(2^{0.34n_{AB}}\right)`$, which is a slight improvement. Let us also notice that this scaling could be further improved by using several levels of nesting, i.e. by replacing the preliminary search over the primary variables by a another nested structured search (see the analysis of the circuit-based counterpart of this idea in the Appendix of J.Cerf et al. (2000)).

## VIII Conclusion

We have introduced a new quantum search algorithm combining the approach based on local adiabatic evolution developed in Roland and Cerf (2002) and the nesting technique introduced in J.Cerf et al. (2000). It allows one to adiabatically solve structured search problems with an improved complexity over a naive adiabatic search that would not exploit the structure of the problem.

The basic idea is to perform a preliminary adiabatic search over a reduced number of variables of the problem in order to keep only a superposition of the assignments that respect the constraints of this partial problem, and then to complete these partial solutions by finding satisfying assignments for the remaining variables. We have seen that to implement this algorithm, the global adiabatic evolution (stage C) has to be discretized, which makes it possible to nest the preliminary search (stages A and B) into the global one. Each step of the algorithm requires to alternate partial adiabatic searches backwards and forwards with global search operations.

A complexity analysis shows that the average computation time of this adiabatic algorithm, although still exponential, grows with a reduced exponent compared to quantum unstructured search algorithms or classical known algorithms to solve a problem such as $`3`$-SAT.

###### Acknowledgements.

J.R. acknowledges support from the Belgian foundation FRIA. This work was funded in part by the Communauté Française de Belgique under grant ARC 00/05-251, by the IUAP programme of the Belgian government under grant V-18, and by the EU under project RESQ (IST-2001-35759).


text:
# THE LAX CONJECTURE IS TRUE

## Abstract

In 1958 Lax conjectured that hyperbolic polynomials in three variables are determinants of linear combinations of three symmetric matrices. This conjecture is equivalent to a recent observation of Helton and Vinnikov.

Key words: hyperbolic polynomial, Lax conjecture, hyperbolicity cone, semidefinite representable
AMS 2000 Subject Classification: 15A45, 90C25, 52A41

Consider a polynomial $`p`$ on $`𝐑^n`$ of degree $`d`$ (the maximum of the degrees of the monomials in the expansion of $`p`$). We call $`p`$ homogeneous if $`p(tw)=t^dp(w)`$ for all real $`t`$ and vectors $`w𝐑^n`$: equivalently, every monomial in the expansion of $`p`$ has degree $`d`$. We denote the set of such polynomials by $`𝐇^n(d)`$. By identifying a polynomial with its vector of coefficients, we can consider $`𝐇^n(d)`$ as a normed vector space of dimension $`\left(\genfrac{}{}{0pt}{}{n+d1}{d}\right)`$.

A polynomial $`p𝐇^n(d)`$ is hyperbolic with respect to a vector $`e𝐑^n`$ if $`p(e)0`$ and, for all vectors $`w𝐑^n`$, the univariate polynomial $`tp(wte)`$ has all real roots. The corresponding hyperbolicity cone is the open convex cone (see )

$$\{w𝐑^n:p(wte)=0t>0\}.$$

For example, the polynomial $`w_1w_2\mathrm{}w_n`$ is hyperbolic with respect to the vector $`(1,1,\mathrm{},1)`$, since the polynomial $`t(w_1t)(w_2t)\mathrm{}(w_nt)`$ has roots $`w_1,w_2,\mathrm{},w_n`$; hence the corresponding hyperbolicity cone is the open positive orthant.

Hyperbolic polynomials and their hyperbolicity cones originally appeared in the partial differential equations literature . They have attracted attention more recently as fundamental objects in modern convex optimization . Three primary reasons drive this interest:

1. the definition of “hyperbolic polynomial” is strikingly simple;
2. the class of hyperbolic polynomials, although not well-understood, is known to be rich — specifically, its interior in $`𝐇^n(d)`$ is nonempty;
3. optimization problems posed over hyperbolicity cones, with linear objective and constraint functions, are amenable to efficient interior point algorithms.

For more details on these reasons, see .

In light of the interest of hyperbolic polynomials to optimization theorists, it is therefore natural to ask: how general is the class of hyperbolicity cones? In particular, do hyperbolicity cones provide a more general model for convex optimization than “semidefinite programming” (the study of optimization problems with linear objectives and constraints and semidefinite matrix variables )?

We begin with some easy observations. A rich source of examples of hyperbolicity cones are semidefinite slices, by which we mean sets of the form

(1)
$$\{w:\underset{j=1}{\overset{n}{}}w_jG_j𝐒_{++}^d\},$$

for matrices $`G_1,G_2,\mathrm{},G_n`$ in the space $`𝐒^d`$ of all $`d`$-by-$`d`$ real symmetric matrices, where $`𝐒_{++}^d`$ denotes the positive definite cone. Such cones are, in particular, “semidefinite representable” in the sense of .

###### Proposition 2

Any nonempty semidefinite slice is a hyperbolicity cone.

Proof Suppose the semidefinite slice (1) contains the vector $`\widehat{w}`$. We claim the polynomial $`p`$ on $`𝐑^n`$ defined by

(3)
$$p(w)=det\underset{j}{}w_jG_j$$

is hyperbolic with respect to $`\widehat{w}`$, with corresponding hyperbolicity described by (1). Clearly $`p`$ is homogeneous of degree $`d`$, and $`p(\widehat{w})>0`$.

Define a matrix $`\widehat{G}=_j\widehat{w}_jG_j𝐒_{++}^d`$, and notice, for any vector $`w𝐑^n`$ and scalar $`t`$, we have

$`p(wt\widehat{w})=det{\displaystyle \underset{j}{}}(w_jt\widehat{w}_j)G_j=det\left({\displaystyle \underset{j}{}}w_jG_jt\widehat{G}\right)`$
$`\text{ }=(det\widehat{G})det\left(\widehat{G}^{1/2}\left[{\displaystyle \underset{j}{}}w_jG_j\right]\widehat{G}^{1/2}tI\right),`$
where $`I`$ denotes the identity matrix. Consequently, the univariate polynomial $`tp(wt\widehat{w})`$ has all real roots, namely the eigenvalues of the symmetric matrix $`H=\widehat{G}^{1/2}[_jw_jG_j]\widehat{G}^{1/2}`$, so $`p`$ is hyperbolic with respect to $`\widehat{w}`$. Furthermore, by definition, $`w`$ lies in the corresponding hyperbolicity cone exactly when these roots (or equivalently, eigenvalues) are all strictly positive. But this property is equivalent to $`H`$ being positive definite, which holds if and only if $`_jw_jG_j`$ is positive definite, as required. $`\mathrm{}`$

The class of semidefinite slices is quite broad. For example, any homogeneous cone (an open convex pointed cone whose automorphism group acts transitively) is a semidefinite slice  (see also ). In particular, therefore, any homogeneous cone is a hyperbolicity cone, a result first observed in .

What about the converse? When is a hyperbolicity cone a semidefinite slice? How general is the class of hyperbolic polynomials of the form (3)?

In considering a general hyperbolic polynomial $`p`$ on $`𝐑^n`$ with respect to a vector $`e`$, we can suppose, after a change of variables, that $`e=(1,0,0,\mathrm{},0)`$ and $`p(e)=1`$. Consider the first nontrivial case, that of $`n=2`$. By assumption, the polynomial $`tp(t,1)`$ has all real roots, which we denote $`g_1,g_2,\mathrm{},g_d`$, so for some nonzero real $`k`$ we have the identity

$$p(t,1)=k\underset{j=1}{\overset{d}{}}(g_jt).$$

By homogeneity, for any vector $`(x,y)𝐑^2`$ with $`y0`$, we deduce

$$p(x,y)=y^dp(\frac{x}{y},1)=y^dk\underset{j=1}{\overset{d}{}}\left(g_j+\frac{x}{y}\right)=k\underset{j=1}{\overset{d}{}}(g_jy+x).$$

By continuity and the fact that $`p(1,0)=1`$, we see

$$p(x,y)=\underset{j=1}{\overset{d}{}}(g_jy+x)=det(xI+yG)$$

for all $`(x,y)𝐑^2`$, where $`G`$ is the diagonal matrix with diagonal entries $`g_1,g_2,\mathrm{},g_d`$. Thus any such hyperbolic polynomial $`p`$ does indeed have the form (3).

What about hyperbolic polynomials in more than two variables? The following conjecture  proposes that all hyperbolic polynomials in three variables are likewise easily described in terms of determinants of symmetric matrices.

###### Conjecture 4 (Lax, 1958)

A polynomial $`p`$ on $`𝐑^3`$ is hyperbolic of degree $`d`$ with respect to the vector $`e=(1,0,0)`$ and satisfies $`p(e)=1`$ if and only if there exist matrices $`B,C𝐒^d`$ such that $`p`$ is given by

(5)
$$p(x,y,z)=det(xI+yB+zC).$$

An obvious consequence of this conjecture would be that, in $`𝐑^3`$, hyperbolicity cones and semidefinite slices comprise identical classes.

A polynomial on $`𝐑^2`$ is a real zero polynomial  if, for all vectors $`(y,z)𝐑^2`$, the univariate polynomial $`tq(ty,tz)`$ has all real roots. Such polynomials are closely related to hyperbolic polynomials via the following elementary result.

###### Proposition 6

If $`p`$ is a hyperbolic polynomial of degree $`d`$ on $`𝐑^3`$ with respect to the vector $`e=(1,0,0)`$, and $`p(e)=1`$, then the polynomial on $`𝐑^2`$ defined by $`q(y,z)=p(1,y,z)`$ is a real zero polynomial of degree no more than $`d`$, and satisfying $`q(0,0)=1`$.

Conversely, if $`q`$ is a real zero polynomial of degree $`d`$ on $`𝐑^2`$ satisfying $`q(0,0)=1`$, then the polynomial on $`𝐑^3`$ defined by

(7)
$$p(x,y,z)=x^dq(\frac{y}{x},\frac{z}{x})(x0)$$

(extended to $`𝐑^3`$ by continuity) is a hyperbolic polynomial of degree $`d`$ on $`𝐑^3`$ with respect to $`e`$, and $`p(e)=1`$.

Proof To prove the first statement, note that for any point $`(y,z)𝐑^2`$ and complex $`\mu `$, if $`q(\mu (y,z))=0`$ then $`\mu 0`$ and $`0=p(1,\mu y,\mu z)=\mu ^dp(\mu ^1,y,z)`$, using the homogeneity of $`p`$. So, by the hyperbolic property, $`\mu ^1`$ is real, and hence so is $`\mu `$. The remaining claims are clear.

For the converse direction, since $`q`$ has degree $`d`$, clearly $`p`$ is well-defined and homogeneous of degree $`d`$ and satisfies $`p(e)=1`$. If $`p(\mu ,y,z)=0`$, then either $`\mu =0`$ or $`q(\mu ^1(y,z))=0`$, in which case $`\mu ^1`$ and hence also $`\mu `$ must be real. $`\mathrm{}`$

(Notice, in the first claim of the proposition, that the polynomial $`q`$ may have degree strictly less than $`d`$: consider, for example, the case $`p(x,y,z)=x^d`$.)

Helton and Vinnikov \[7, p. 10\] observe the following result, based heavily on .

###### Theorem 8

A polynomial $`q`$ on $`𝐑^2`$ is a real zero polynomial of degree $`d`$ and satisfies $`q(0,0)=1`$ if and only if there exist matrices $`B,C𝐒^d`$ such that $`q`$ is given by

(9)
$$q(y,z)=det(I+yB+zC).$$

(Notice, as in the Lax conjecture, the “if” direction is immediate.)

We claim that Theorem 8 is equivalent to the Lax conjecture. To see this, suppose $`p`$ is a hyperbolic polynomial of degree $`d`$ on $`𝐑^3`$ with respect to the vector $`e=(1,0,0)`$, and $`p(e)=1`$. Then by Proposition 6, the polynomial on $`𝐑^2`$ defined by $`q(y,z)=p(1,y,z)`$ is a real zero polynomial of degree $`d^{}d`$, and satisfying $`q(0,0)=1`$. Hence by Theorem 8, equation (9) holds: we can assume $`d^{}=d`$ by replacing $`B,C𝐒^d^{}`$ with block diagonal matrices $`\text{Diag}(B,0),\text{Diag}(C,0)𝐒^d`$. Then, by homogeneity, for $`x0`$,

$`p(x,y,z)=x^dp(1,{\displaystyle \frac{y}{x}},{\displaystyle \frac{z}{x}})=x^dq({\displaystyle \frac{y}{x}},{\displaystyle \frac{z}{x}})`$
$`=x^ddet\left(I+{\displaystyle \frac{y}{x}}B+{\displaystyle \frac{z}{x}}C\right)=det(xI+yB+zC).`$
as required. The converse direction in the Lax conjecture is immediate.

Conversely, let us assume the Lax conjecture, and suppose $`q`$ is a real zero polynomial of degree $`d`$ on $`𝐑^2`$ satisfying $`q(0,0)=1`$. (The converse direction in Theorem 8 is immediate.) Then by Proposition 6 the polynomial $`p`$ defined by equation (7) is a hyperbolic polynomial of degree $`d`$ on $`𝐑^3`$ with respect to $`e`$, and $`p(e)=1`$. According to the Lax conjecture, equation (5) holds, so

$$q(y,z)=p(1,y,z)=det(I+yB+zC),$$

as required. $`\mathrm{}`$

The exact analogue of the Lax conjecture fails in general for polynomials in $`n>3`$ variables. To see this, note that the set of polynomials on $`𝐑^n`$ of the form $`wdet_jw_jG_j`$ (where $`G_1,G_2,\mathrm{},G_n𝐒^d`$) has dimension at most $`n\left(\genfrac{}{}{0pt}{}{d+1}{2}\right)`$, being an algebraic image of a vector space of this dimension. If the degree $`d`$ is large, this dimension is certainly smaller than the dimension of the set of hyperbolic polynomials: as we observed above, this latter set has nonempty interior in the space $`𝐇^n(d)`$ (by a result of Nuij \[6, Thm 2.1\]), and so has dimension $`\left(\genfrac{}{}{0pt}{}{n+d1}{d}\right)`$.

More concretely, consider the polynomial defined by $`p(w)=w_1^2_2^nw_j^2`$ for $`w𝐑^n`$. This polynomial is hyperbolic of degree $`d=2`$ with respect to the vector $`(1,0,0,\mathrm{},0)`$, and yet cannot be written in the form $`det_jw_jG_j`$ for matrices $`G_1,G_2,\mathrm{},G_n𝐒^2`$ if $`n>3`$. To see this, choose any nonzero vector $`w`$ satisfying $`w_1=0`$, and such that the first row of the matrix $`_jw_jG_j`$ is zero.

The question of whether all hyperbolicity cones are semidefinite slices, or, more generally, are semidefinite representable, appears open.

Acknowledgement We are very grateful to the Institute for Mathematics and its Applications at the University of Minnesota for their hospitality during our work on this topic.