text:
# A model for fluctuating inflaton coupling: (s)neutrino induced adiabatic perturbations and non-thermal leptogenesis

\[

## Abstract

We discuss an unique possibility of generating adiabatic density perturbations and leptogenesis from the spatial fluctuations of the inflaton decay rate. The key assumption is that the initial isocurvature perturbations are created in the right handed sneutrino sector during inflation which is then converted into adiabatic perturbations when the inflaton decays. We discuss distinct imprints on the cosmic micro wave background radiation, which can distinguish non-thermal versus thermal leptogenesis.

\]

Inflation is the main contender for explaining the observed adiabatic density perturbations with a nearly scale invariant spectrum . However, recently various alternative mechanisms for generating adiabatic density perturbations have been discussed, particularly converting the isocurvature perturbations of some light field into the adiabatic perturbations in the post-inflationary Universe . Another interesting proposal is that the perturbations could be generated from the fluctuations of the inflaton coupling to the Standard Model degrees of freedom . It has been argued that the inflaton coupling strength to the ordinary matter, instead of being a constant, could depend on the vacuum expectation values (VEV) of various fields in the theory. These fields are none other than the flat directions of the minimal supersymmetric standard model (MSSM). However the authors in  treated the flat directions without considering the fact that the flat directions are lifted at a non-renormalizable level. In , the authors have demonstrated the importance of non-renormalizable potential terms for the flat directions, which leads to dramatic changes in the estimation of the amplitude of the density perturbations in the original scheme. The amplitude of the perturbations dampens after the end of inflation, because the flat directions evolve after the end of inflation until the decay of the inflaton. The damping of the amplitude of the perturbations acts as a main challenge for realizing such a novel scheme (for a review on MSSM flat direction and cosmology, see ).

The idea is that if the MSSM condensates (made up of squarks and sleptons) are light during inflation then their quantum fluctuations can give rise to spatial fluctuations in the inflaton coupling strength. When the inflaton decays, the adiabatic density perturbations are created because the isocurvature perturbations generated by the flat direction during inflation is transferred to the adiabatic ones right at the time of decay. If the flat direction evaporates into baryons, it will give rise to the baryon isocurvature fluctuations, which can be constrained from cosmic microwave background radiation (CMB).

A particularly interesting implementation of this scenario is to consider the right handed (s)neutrinos (supersymmetry guarantees the presence of the right handed sneutrino). If the right handed (s)neutrino is Majorana, then it provides a natural explaining for the observed small neutrino masses via see-saw mechanism, $`m_\nu \left(m_\mathrm{D}^2/M_N\right)`$ , where $`m_\mathrm{D}`$ is the Dirac mass obtained from the Higgs VEV.

The right handed (s)neutrino is also a source for $`L`$ and/or $`BL`$ violation. Therefore it can be responsible for the observed baryon asymmetry. Leptogenesis requires $`L`$ or $`BL`$ violating interactions, $`C`$ and $`CP`$ asymmetry, and an out of equilibrium condition. The first two are well served by the right handed (s)neutrinos, and the last condition naturally arises in any inflationary cosmology.

Now let us briefly discuss the status of inflation. Undoubtfully, inflation is the most natural mechanism which makes the Universe homogeneous, flat, and isotropic. A single field slow roll inflation is also the most beautiful way of explaining adiabatic density perturbations. However until now there is hardly any connection between the inflationary sector and the known particle physics sector. In most of the cases the inflaton is a gauge singlet, which leads to some degree of speculation on how the inflaton couples to the SM gauge particles ? what is the inflaton potential ?, etc. The coupling of the inflaton to the SM fields is essential if the inflationary paradigm wishes to make connection to the hot big bang cosmology.

In this paper our aim is to present a simple toy model where we illustrate two important aspects. The first one is to consider that the sneutrino induced isocurvature fluctuations which can generate adiabatic perturbations through inflaton decay. Second, we will show that in a non-thermal leptogenesis scenario the produced baryon-isocurvature fluctuations can be testable from CMB.

In this regard it is natural to come up with a scenario where the inflaton couples only through the right handed (s)neutrino sector. Note that the right handed (s)neutrino acts as a mediator which connects the two disparate sectors, e.g., the inflaton and the SM via lepton (sleptons) and Higgs (Higgsinos). Therefore such a model is not only economical in terms of achieving density perturbations and lepton asymmetry, but also providing us with the correct relativistic species required for the nucleosynthesis.

For the purpose of illustration, let us consider a very simple toy model,

$$W\frac{1}{2}g\mathrm{\Phi }\mathrm{𝐍𝐍}+h\mathrm{𝐍𝐇}_u𝐋+\frac{1}{2}M_N\mathrm{𝐍𝐍},$$
(1)
where $`\mathrm{\Phi }`$, $`𝐍`$, $`𝐋`$, and $`𝐇`$, respectively stand for the inflaton, the right handed neutrino, the lepton doublet, and the Higgs (which gives mass to the top quark) superfields. Also $`M_N`$ denotes right handed (s)neutrino masses, and $`g,h`$ correspond to the Yukawas. Being a SM gauge singlet the inflaton can naturally couple to the right handed neutrino sector with a renormalizable and a non-renormalizable coupling

$$g=g_0\left(1+\frac{𝐍}{M_\mathrm{p}}+\mathrm{}\right),$$
(2)
where we assume that the non-renormalizable scale is the Planck mass, $`M_\mathrm{p}=2.4\times 10^{18}`$ GeV. For simplicity we have omitted all the indices in $`h`$ matrix and superfields, and we work on a basis where the Majorana mass matrix is diagonal. Further simplifications can be made for almost degenerate right handed (s)neutrinos where $`M_N`$ is essentially the same for all. It is also conceivable in this case that the inflaton is coupled with the same strength to three right handed (s)neutrinos with a mass hierarchy $`M_Nm_\varphi `$, where $`m_\varphi `$ can be treated as the mass of the inflaton in the minimum of its potential. Note that this is just a working example of non-thermal leptogenesis. We will highlight why we stress upon non-thermal leptogenesis compared to thermal leptogenesis. More complicated scenarios on non-thermal leptogenesis can be constructed .

The inflaton sector is still unknown, except that it is responsible for driving inflation, which could be for e.g., brane driven inflation, fast rolling inflation, kinetic driven inflation, assisted inflation, false vacuum inflation, etc. We further assume that the inflaton decays perturbatively. At this point one might suspect that the above coupling, $`g`$, in Eq. (1), would give rise to a large mass contribution to the sneutrino through the inflaton vev, and the sneutrino would simply roll down to the bottom of the potential. However note that the inflaton vev need not be always large in order to inflate, for instance in a false vacuum inflation the inflaton could be fairly close to the bottom of its own potential. Nevertheless, let us imagine that we are working in a regime where the coupling, $`g`$, is such that the effective mass for the sneutrino is less than the Hubble expansion, $`m_{eff,\varphi }<M_{eff,N}H`$, where $`m_{eff,\varphi }`$ is the effective mass of inflaton during inflation. We will comment on the situation when the sneutrino masses are heavy compared to the Hubble expansion during inflation. From here onwards we remove the subscript, $`eff`$.

An important point to note here is that during inflation the quantum fluctuations are created in the sneutrino sector, and the perturbations in the inflaton sector is assumed to be negligible, therefore, the perturbations arises only from the known particle physics sneutrino sector. The perturbations on a comoving scale larger than the Hubble scale can be foliated in terms of the curvature perturbations on a finite energy density surface: $`ds^2=a^2(t)\left(1+2\zeta \right)dx^idx^j`$, where $`\zeta `$ is a metric perturbation written in a proper coordinate system. In presence of more than one scalar fields the total curvature perturbations $`\zeta `$ evolves outside the horizon due to non-vanishing pressure perturbations .

$$\dot{\zeta }=\frac{H}{\rho +P}\delta P.$$
(3)
where $`\rho ,P`$ are the energy density and the pressure. For a single field inflation, $`\zeta =\mathrm{constant}`$, but in a multi-field case $`\delta P`$ is a non-zero quantity due to the entropy perturbations, which can be defined in our case as

$$S_{\varphi ,\stackrel{~}{N}}=3\left(\zeta _\varphi \zeta _{\stackrel{~}{N}}\right)=3H\left(\frac{\delta \rho _\varphi }{\dot{\rho }_\varphi }\frac{\delta \rho _{\stackrel{~}{N}}}{\dot{\rho }_{\stackrel{~}{N}}}\right).$$
(4)
Over-dot denotes differentiation w.r.t. coordinate time. Following our assumption the initial entropy perturbation becomes $`S_{\varphi ,\stackrel{~}{N}}3\zeta _{\stackrel{~}{N}}`$. For the Gaussian perturbations, we obtain

$$𝒫_{\stackrel{~}{N}}^{1/2}=\frac{H_{}}{2\pi },$$
(5)
where $``$ denotes when the interesting perturbations leave the horizon, $`k=a_{}H_{}`$. Note that the entropy perturbations feed the total curvature perturbations, therefore the entropy perturbations along with the individual perturbations, $`\zeta _\varphi ,\zeta _{\stackrel{~}{N}}`$, evolve in time. Though we do not prove this here, but intuitively we can see that in order to obtain the adiabatic density perturbations, the total curvature perturbation must become constant outside the horizon at the time of inflaton decay, when $`\zeta _\varphi |_{decay}\zeta _{\stackrel{~}{N}}`$, and therefore, $`\zeta `$, becomes constant on large scales. Thus converting its initial isocurvature fluctuations into the adiabatic ones.

Note that besides the fluctuations in the sneutrino sector, there will be fluctuations in the Higgses and the sleptons also. However their perturbations will not account for the baryon isocurvature fluctuations in the above set up, for the time being we will neglect them.

The spectral index for the perturbations can be written as

$$n_\zeta 1\frac{d\mathrm{ln}𝒫_{\stackrel{~}{N}}}{d\mathrm{ln}k}=2\frac{\dot{H}_{}}{H_{}^2}+\frac{2}{3}\frac{M_N^2}{H_{}^2}.$$
(6)
Therefore as long as the Hubble expansion is slowly varying, and $`M_NH`$, we can obtain a scale invariant density perturbation.

Now let us study the decay of the inflaton. The main decay mode of the inflaton is four-body final states consisting of two Higgs/Higgsino-lepton/slepton particles (and their $`CP`$ transforms), this is due to the fact that the inflaton is decaying via off-shell (s)neutrino. The effective superpotential after integrating out $`𝐍`$, is given by

$$W_{eff}\frac{1}{2M_N^2}gh^2\mathrm{\Phi }(𝐇_u𝐋)(𝐇_u𝐋).$$
(7)
For simplicity we may consider a situation when the (s)neutrinos are almost degenerate, e.g. $`\mathrm{\Delta }M_N<M_N`$, and the Yukawa texture is such that the diagonal entries (h) and off-diagonal entries ($`h^{}`$) follow $`h^{}<h`$.

There are total nine final states, seven of them consist of two fermions and two scalars, and there are also two final states consisting of four scalars. Summing up all the final states, the decay rate and the final reheat temperature are given by

$$\mathrm{\Gamma }_\mathrm{d}\frac{21g^2h^4m_\varphi ^5}{2^{14}\pi ^5M_N^4},\frac{T_\mathrm{R}}{m_\varphi }\frac{10^{7/2}gh^2m_\varphi ^{3/2}M_\mathrm{P}^{1/2}}{M_N^2}.$$
(8)
However note that the reheat temperature obtains a spatial fluctuations due to Eq. (2),

$$\frac{\delta T_\mathrm{R}}{T_\mathrm{R}}=\frac{1}{3}\frac{\delta g}{g}\frac{\delta \stackrel{~}{N}}{3M_\mathrm{p}}\frac{H_{}}{6\pi M_\mathrm{p}}.$$
(9)
The factor $`1/3`$ arises because during the decay of the inflaton the average energy density goes as $`\rho a^3`$, see for details . For the Gaussian perturbations, $`\stackrel{~}{N}H`$, and following Eq. (5) we obtain the right amplitude for the density perturbations provided $`H_{}10^5M_\mathrm{p}`$. Note that there is no damping in the sneutrino fluctuations after the end of inflation in this case. Further note that if $`\stackrel{~}{N}H`$, then the Gaussian amplitude of the perturbations will be damped , which we would like to avoid in this example.

Let us now obtain the lepton asymmetry in this model. The $`CP`$ asymmetry is obtained through the inference between the tree level and the one loop (vertex and self energy) corrected diagrams, Net $`CP`$ asymmetry in the off-shell case is quite different compared to the on-shell leptogenesis . The self energy correction comes out to be twice as much as the vertex correction for $`m_\varphi M_N`$. Final $`CP`$ asymmetry is then given by

$$ϵ_{CP}\frac{3}{8\pi }\times \frac{\underset{i,n,l}{}\frac{\mathrm{Im}\left[(\mathrm{𝐡𝐡}^{})_{ni}(\mathrm{𝐡𝐡}^{})_{nl}(\mathrm{𝐡𝐡}^{})_{il}\right]m_\varphi ^2}{M_i^3M_n^2M_l}}{_{i,n}\frac{([\mathrm{𝐡𝐡}^{}]_{in})^2}{M_i^2M_n^2}},$$
(10)
where $`i,n,l=1,2,3`$. The produced lepton asymmetry reads as

$$\frac{n_L}{n_\varphi }\frac{3}{\pi }\frac{\delta hh^2}{h}\frac{\mathrm{\Delta }M_N}{M_N}\left(\frac{m_\varphi }{M_N}\right)^2,$$
(11)
where $`\delta h`$ is nearly equal diagonal entries of the Yukawa matrix. For nearly degenerate case $`\delta h/h\mathrm{\Delta }M_N/2M_N`$.

The total asymmetry in the baryons (after taking into account of the sphaleron effects) can be expressed as

$`\eta _\mathrm{B}`$ $`=`$ $`\left({\displaystyle \frac{n_B}{n_\varphi }}\right)\left({\displaystyle \frac{n_\varphi }{s}}\right)`$ (12)
$``$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle \frac{\delta hh^2}{h^3}}{\displaystyle \frac{\mathrm{\Delta }M_N}{M_N}}\left({\displaystyle \frac{M_Nm_\nu }{H_u^0^2}}\right)\times \left({\displaystyle \frac{m_\varphi }{M_N}}\right)^2\left({\displaystyle \frac{T_\mathrm{R}}{m_\varphi }}\right),`$ (13)
where $`s=(2\pi ^2/45)g_{}T_\mathrm{R}^3`$. Here $`n_\varphi /s`$ denotes the dilution from reheating. By using the expression for the reheat temperature and $`m_\nu (h^2H_u^0^2/M_N)`$, we finally obtain

$$\eta _\mathrm{B}4.10^{49/2}g\frac{\delta hh^2}{h^3}\frac{\mathrm{\Delta }M_N}{M_N}\frac{m_\varphi ^{7/2}M_\mathrm{P}^{1/2}}{M_N^2H_u^0^4}(1\mathrm{G}\mathrm{e}\mathrm{V})^2,$$
(14)
where we have considered $`m_\nu 0.1`$ eV, and $`H_u^0=174`$ GeV. If we demand degenerate light neutrino masses, then we can further simplify, Eq. (14), $`\eta _\mathrm{B}`$

$$\frac{2.10^{49/2}gh^2}{h^2}\left(\frac{\mathrm{\Delta }M_N}{M_N}\right)^2\frac{m_\varphi ^{7/2}M_\mathrm{P}^{1/2}}{M_N^2H_u^0^4}(1\mathrm{G}\mathrm{e}\mathrm{V})^2.$$
(15)
Some numerical examples for nearly degenerate heavy right handed (s)neutrinos, with $`M_N=10m_\varphi `$ and $`10^1h^{}/h1`$, we obtain the desired baryon asymmetry for $`10^3g1`$ and $`10^{11}\mathrm{GeV}m_\varphi 10^{13}`$ GeV, which result in reheat temperature: $`10^6\mathrm{GeV}T_\mathrm{R}10^8`$ GeV. At this point one might worry upon the coupling strength, $`g`$, because the inflaton picks up an effective mass term, $`g\stackrel{~}{N}`$, which has to be smaller than the inflaton mass, $`m_\varphi `$, arising solely from the inflaton sector, in order to keep the successes of slow roll inflation. On the other hand in order to generate the Gaussian fluctuations, $`\stackrel{~}{N}H_{}10^5M_p`$, see Eq. (9). This leads to $`g1`$, however, an exact magnitude of $`g`$ will depend on a particular inflationary model. Further note that the reheat temperature is well below thermal and non-thermal gravitino over-production .

The most important point is to note that the baryon asymmetry is proportional to $`g`$, see our final result, Eq. (15). Therefore baryons also feel the spatial fluctuations.

$$\frac{\delta \eta _B}{\eta _B}\frac{1}{3}\frac{\delta g}{g}\frac{\delta \stackrel{~}{N}}{3M_\mathrm{p}}\frac{\delta T_\mathrm{R}}{T_\mathrm{R}}0.$$
(16)
The origin of $`1/3`$ factor has the same origin as in Eq. (9). Note that the fluctuations in the baryon asymmetry is proportional to the fluctuations in the inflaton coupling, and therefore fluctuations in the reheat temperature. This shows that the baryonic asymmetry does not follow honest to God the adiabatic density perturbations, instead perturbation in baryons is correlated baryon-isocurvature in nature.

The baryon-isocurvature fluctuation leaves its imprint on cosmic micro wave background radiation. Moreover the fluctuations are correlated. The WMAP data provides mild constraint on the correlated-cold dark matter-isocurvature fluctuations , which can be translated in terms of the baryon isocurvature fluctuations as $`|S_B/\zeta |<0.32(\mathrm{\Omega }_{CDM}/\mathrm{\Omega }_B)1.85`$ at $`95\%`$ confidence level, where $`S_B`$ is the baryon-isocurvature fluctuations. In our toy model $`S_B=\delta \eta _B/\eta _B=\delta T_\mathrm{R}/T_\mathrm{R}`$. In particular $`\zeta =H\delta \rho /\dot{\rho }=(1/4)\delta \rho _\gamma /\rho _\gamma =\delta T_\mathrm{R}/T_\mathrm{R}`$, where the subscript $`\gamma `$ denotes radiation. Therefore we find $`|S_B/\zeta |=1`$ in our case, which is well within the WMAP constraint on the baryon-isocurvature perturbations.

We point out here that the above feature of correlated baryon-isocurvature perturbations is only present in a non-thermal case, this is the reason why we pursued on non-thermal leptogenesis. In non-thermal leptogenesis there is an explicit dependence on the reheat temperature, see Eq. (12). This is indeed an interesting feature of a non-thermal leptogenesis which is absent in a thermal case. In a thermal leptogenesis the net asymmetry is proportional to a $`CP`$ asymmetry and not to a temperature 111In general thermal/non-thermal leptogenesis provides the net asymmetry as $`\eta _{B/L}ϵ_{CP}\times f(T_d/M)`$, where $`T_D`$ is the temperature of the decaying particle, e.g. the right handed neutrino. In case of thermal leptogenesis $`T_dM`$, and there is no temperature dependence. However in a non-thermal case $`T_dM`$.. We find that the constraints on the baryon-isocurvature perturbations can act as a tool for differentiating thermal versus non-thermal leptogenesis mechanisms.

Before we conclude our paper, we comment on couple of interesting points. In an opposite limit when $`m_\varphi >M_N`$, the inflaton decays via on-shell right handed (s)neutrino to the SM leptons and Higgs. This case is even better because during inflation the condition $`\stackrel{~}{N}H_{}`$ is satisfied even better, because the sneutrino is lighter, $`M_Nm_\varphi H`$. However one has to ensure that a thermal regeneration of the baryon asymmetry is really small. Finally we comment on a heavy right handed neutrino masses compared to the Hubble expansion, in this case the sneutrino perturbations will be $`\chi ^2`$ in nature, and usually the amplitude of the perturbations comes out to be small . Nevertheless, they can also provide interesting imprints on CMB through the tilt in the spectral index. We leave this for future investigation.

In fact we could also imagine perturbing the other Yukawa coupling, $`h`$, similar to Eq. (2), and the inflaton mass due to the fluctuations in the sneutrino vev. By inspecting the reheat temperature, Eq. (8), and the baryon asymmetry, Eq. (12), we obtain $`\delta T_\mathrm{R}/T_\mathrm{R}=(2/3)(\delta h/h)=\delta \eta _\mathrm{B}/\eta _\mathrm{B}`$, assuming that the fluctuations are arising only from the diagonal elements of the Yukawas. In this case the prediction on the baryon-isocurvature fluctuations remains, $`|S_\mathrm{B}/\zeta |=1`$. However the fluctuating inflaton mass gives rise to $`|S_\mathrm{B}/\zeta |=1.4`$, which is still within the WMAP limit on $`|S_\mathrm{B}/\zeta |<1.85`$ @ $`95\%`$ confidence level .

As a final remark, there could be other sources for the isocurvature perturbations during inflation, including the most competitive candidate “cold dark matter”. However within supersymmetry excellent conditions arise naturally for their thermal production.

To summarize our paper, we point out that any supersymmetric leptogenesis scenario is a potential candidate for succeeding in generating the adiabatic density perturbations from the sneutrino fluctuations during inflation and generating the baryon asymmetry. The nature of perturbations (Gaussian/non-Gaussian) certainly depends on the mass scales, e.g., for $`M_NH_{}`$, and $`\stackrel{~}{N}H_{}`$ the perturbations are Gaussian. This may not be the case if $`M_N=3H`$ during inflation, see . Note that our model is economical because it achieves several goals at a time. Finally we have found a very important bench mark which can potentially differentiate thermal versus non-thermal leptogenesis from CMB.

A. M. is a CITA National fellow. The author has benefited from the discussion with Rouzbeh Allahverdi, Robert Brandenberger, Zurab Berezhiani, Guy Moore, David Wands, and valuable comments from the anonymous referees.


format: text

text:
# The Evolution of Protoplanetary Disk Edges

## 1 Introduction

The discovery of extra-solar planetary systems (Mayor & Queloz, 1995) has driven a flood of renewed interest in planetary formation and disk evolution models. The focus of many models has been on mechanisms that couple the planets to the evolution of protostellar disks, the tidal interactions between a planet and a protoplanetary disk (Artymowicz & Lubow, 1994; Ogilvie & Lubow, 2002) and resulting planetary migration (Lin & Papaloizou, 1986; Trilling et al., 1998; Nelson et al., 2000; Murray et al., 1998). Gaps and holes in protostellar disks in some cases have been detected through the study of spectral energy distributions (Rice et al., 2003; Calvet, 2002; Koerner, Sargent & Beckwith, 1993) and direct imaging (Jayawardhana et al., 1998; Augereau et al., 1999; Weinberger et al., 1999; Grady et al., 2001). Since structure in protostellar disks can arise through the gravitational interaction between the disk and planets, detailed studies of these disks may allow us to place constraints on the properties of the planets in them and on the evolution of these systems.

The formation of gaps in gaseous or planetesimal disks has a long history of theoretical study. It has long been postulated that an external time dependent gravitational potential, such as from a satellite or planet, could resonantly drive waves at Lindblad resonances (LRs) in a gaseous or planetesimal disk (Goldreich & Tremaine, 1978, 1979; Lin & Papaloizou, 1979, 1993; Ward, 1997). Because these waves carry angular momentum, they govern the way that planets open gaps in disks (e.g., Bryden et al. 1999; Artymowicz & Lubow 1994) and can result in the orbital migration of planets (e.g., Nelson et al. 2000). Previous works have derived and numerically verified conditions for a planet to open a gap in a gas disk (e.g., Bryden et al. 1999). These studies focused on when gaps occur and explored the accretion rate through the gap as material was cleared to either side. Emphasis was also placed on the relationship between gap properties, accretion and orbital migration rates (e.g., Kley 1999; Nelson et al. 2000).

However, little numerical work has focused on the expected sizes and depths of the gaps in different physical regimes. Such studies are warranted because, as discussed above, gaps and holes in protoplanetary disks are detectable and so will in principle allow us to place constraints on the properties of the planets that are responsible for maintaining them.

The physical properties of protoplanetary disks should vary as the disk ages. We parametrize the disk viscosity, $`\nu `$, with a Reynolds number, $`r^2\mathrm{\Omega }/\nu `$ where $`\mathrm{\Omega }`$ is the angular rotation rate of a particle in a circular orbit around the star and $`r`$ is the radius. Using an $`\alpha `$ form for the viscosity,

$$=\alpha ^1\left(\frac{v_c}{c_s}\right)^2=\alpha ^1\left(\frac{r}{H}\right)^2$$
(1)
where $`c_s`$ is the sound speed, $`H`$ is the vertical scale height of the disk and $`v_c=r\mathrm{\Omega }`$ is the velocity of a particle in a circular orbit. We expect that the sound speed and vertical scale height of a cooling gaseous disk will drop with time, increasing the Reynolds number. For a planetesimal disk, the effective viscosity depends upon the opacity, $`\tau `$, (or area filling factor) of the disk with $`\alpha \frac{\tau }{2}`$ in the limit of low $`\tau `$ (e.g., Murray & Dermott 1999 and citations within). As planetesimals coagulate or are scattered out of the disk, we expect a decrease in the viscosity and so an increase in the Reynolds number.

When viscous forces are overcome by the angular momentum flux driven by density waves, a gap in the disk near the planet opens. So as protoplanetary disks age (and their viscosity drops) we expect that gaps caused by planets will widen and accretion onto these planets will cease (Bryden et al., 1999; Kley, 1999). Because of the longer timescale associated with the end of accretion, it may be that forthcoming surveys will be dominated by partially cleared disks. The advent of new instruments such as SIRTF, ALMA, SMA and VLTI make a study of this limit particularly interesting and relevant (e.g., Wolf et al. 2002).

In this paper we numerically investigate the properties of gaps, exploring a larger range of parameter space than have previous studies. In particular we extend our simulations toward the low viscosity limit, which may govern a longer timescale in the evolution of protoplanetary disks. We critically review our results in the context of previous theories of gap formation and attempt to extend these models to embrace our new results. In a follow up paper we will explore the observational properties of our simulated disks (Varnière et al., 2004).

The structure of our paper is as follows. In section 2 we review the theory of gap formation. In section 3 we present our numerical methods and the input conditions for our models. In section 4 we present our results and in section 5 we present a critical discussion of the results along with attempts to extent previous gap formation models. A discussion and conclusion follows in section 6.

## 2 Previous Theory

Linear theory has been developed to predict the angular momentum flux driven by spiral densities waves exited by a planet, (Goldreich & Tremaine, 1978; Donner, 1979; Lin & Papaloizou, 1979; Artymowicz, 1993). This theory has successfully been used to predict the formation of gaps seen in numerical simulations (Lin & Papaloizou, 1993; Bryden et al., 1999). If a gap is opened then the angular momentum flux resulting from spiral density waves can be balanced by the angular momentum flux caused by the viscosity of the disk. To theoretically account for these phenomena, we require expressions for the angular momentum flux due to both processes.

As have previous works, we assume that the perturbing gravitational potential from the planet can be expanded in Fourier components

$$\mathrm{\Phi }(r,\theta ,t)=\underset{m}{}W_m(r)\mathrm{exp}\left(im(\theta \mathrm{\Omega }_pt)\right)$$
(2)
where $`\mathrm{\Omega }_p`$ is the angular rotation rate of the pattern, and $`r,\theta `$ are polar coordinates in the plane containing the planet. Here $`\mathrm{\Omega }_p`$ is the angular rotation rate of a planet which we assume is in a circular orbit. The torque $`T_m`$ exerted at the $`m`$-th Lindblad resonance into a disk of surface density $`\mathrm{\Sigma }`$ was found to depend solely on the radial profile of the perturbing potential and be independent of the process of dissipation providing that the waves are completely dissipated (Goldreich & Tremaine, 1978; Donner, 1979; Lin & Papaloizou, 1979; Artymowicz, 1993),

$$T_m=m\pi ^2\mathrm{\Sigma }\frac{|\mathrm{\Psi }_{GT}|^2}{rdD_m/dr}$$
(3)
where

$$\mathrm{\Psi }_{GT}=r\frac{dW_m}{dr}+\frac{2m\mathrm{\Omega }}{m(\mathrm{\Omega }\mathrm{\Omega }_p)}W_m$$
(4)
and

$$D_m=\kappa ^2m^2(\mathrm{\Omega }\mathrm{\Omega }_p)^2.$$
(5)
Here $`\mathrm{\Omega }`$ and $`\kappa `$ are the angular rotation rate and epicyclic frequency of an unperturbed particle in a circular orbit and are functions of radius. The above torque results from $`m`$-armed waves driven at the $`m`$th Lindblad resonance which is located at a radius determined by $`D_m(r_{LR})=0`$. Equation(3) is evaluated at this radius.

For a planet in a circular orbit, the Fourier components of the potential exterior to the planet are

$$W_m=\frac{GM_p}{r}b_{1/2}^m(r_p/r)$$
(6)
where $`b_{1/2}^m`$ is the Laplace coefficient (e.g., see Murray & Dermott 1999) and $`r_p`$ is the radius of the planet from the central star. We write equation (3) in the form

$$T_m=f_m\mathrm{\Sigma }r_p^4\mathrm{\Omega }_p^2\left(\frac{q}{10^3}\right)^2$$
(7)
where $`f_m`$ are unitless constants, $`qM_p/M_{}`$ is the mass ratio of the planet to the central star, and $`\mathrm{\Sigma }`$ is evaluated at the Lindblad resonance. Evaluating equation(3) using summations for the Laplace coefficients (e.g., Murray & Dermott 1999) we find that $`f_1=10^{5.1}`$, $`f_2=10^{4.3}`$, $`f_3=10^{3.7}`$ and $`f_4=10^{3.4}`$. These values will be used in section 4 in our discussion on clearing timescales and in section 5 in our discussion on the edge profile.

If the planet is sufficiently massive and the disk viscosity sufficiently low, the torque caused by the excitation of spiral density waves will overcome the inflow due to viscous accretion in the disk and a gap will open. The angular momentum flux or torque transferred through a radius $`r`$ in a Keplerian viscous disk with constant viscosity and density is given by

$$T_\nu =3\pi \mathrm{\Sigma }\nu r^2\mathrm{\Omega }.$$
(8)
Previous estimates of gap opening criteria have applied equation(3) in the limit for large $`m`$ so that the resonances closest to the planet are taken into account (Lin & Papaloizou, 1979; Bryden et al., 1999). In these cases the total torque resulting from summing the resonant terms from each outer Lindblad resonance (OLR) within a distance $`\mathrm{\Delta }`$ from the planet yields

$$T0.23q^2\mathrm{\Sigma }r_p^2\mathrm{\Omega }^2\left(\frac{r_p}{\mathrm{\Delta }}\right)^3$$
(9)
where $`r_p`$ is the radius of the planet (Lin & Papaloizou, 1993). In what follows we take $`\mathrm{\Delta }`$ to be the distance from the planet to the gap outer boundary. The expression above is valid when the gap is larger than the maximum value of the disk scale height ($`\mathrm{\Delta }>H`$). When the torque from the dissipation of the outward going density waves is balanced by the torque from viscous accretion (equation 8) one can derive a criterion for gap formation, i.e. $`q>40/`$ (Bryden et al., 1999).

Once a gap is opened, the torque from viscous accretion must be balanced by the total torque from the spiral density waves driven at Lindblad Resonances exterior to the gap. As the mass of the planet increases and the viscosity decreases, the gap edges moves away from the planet so that the total torque contains fewer resonant terms. Past the $`m=2`$ OLR at $`r1.3r_p`$ only the weaker $`m=1`$ OLR at $`r1.6r_p`$ can drive spiral density waves. Previous works have suggested that the location of the $`m=2`$ OLR resonance would set the radius of disk edges in low viscosity disks harboring large mass planets (Lin & Papaloizou, 1993). It is possible that as the gap is cleared, gas could accumulate near this resonance (Lin & Papaloizou, 1993). Alternatively, spiral density waves could be driven at a LR via a ‘virtual effect’ as long as the resonance was within a particular distance of the disk edge. This argument has been used to provide a possible explanation for $`m=2`$ structure seen in simulations of binaries (e.g., Lin & Papaloizou 1993, 1979; Savonije et al. 1994). A large mass ratio, such as would be found in a stellar binary system, or extremely small disk viscosity, is required to open a gap out to the $`m=1`$ OLR in a circumbinary disk (Artymowicz & Lubow, 1994).

The balance between torque from spiral density waves and $`T_\nu `$ should allow one to derive a scaling relation for the gap width. If we naively use equation (9) and assume a Keplerian rotation profile, we find

$$\mathrm{\Delta }=0.29q^{2/3}^{1/3}r_p.$$
(10)
While previous numerical studies have confirmed the condition for opening a gap, a general scaling relation to predict the gap width has not yet been developed. However, it is possible to extract some information from published data in order to probe the scaling of gap width with planet mass and disk viscosity. In Figure 1 we present a plot of equation(10) along with data from Bryden et al. (1999) as well as data from our own simulations. The points by Bryden et al. (1999) were taken for their $`q=10^3`$ simulations and the gap widths were measured from their density profile plots. The gap width was taken to be the distance from the planets orbital radius to the location along the outer gap wall where the density is half of that outside the gap. Reynolds numbers were computed using their quoted values of the viscosity parameter $`\alpha _{ss}`$ and

$$\nu =\alpha _{ss}(H/r)^2r^2\mathrm{\Omega }.$$
(11)
For the zero viscosity case of Bryden et al. (1999), we use their quoted value of the numerical viscosity yielding $`\alpha _{ss}=10^5`$.

Inspection of Figure 1 shows that the scaling relation given in equation (10) produces the correct order of magnitude gap width. Equation(10), since it is based upon a high $`m`$ integration, is correct for lower planet masses and higher disk viscosities, providing an explanation for the correct positions of the points in Figure 1 on the high viscosity end. Given that the gap edge is not sharp, there is some ambiguity in the definition of $`\mathrm{\Delta }`$ and hence in the exact comparison between equation(10) and the simulations. We take up this point again in section 4. In spite of this uncertainty, the failure of the scaling relation to predict the gap width across a wide range of Reynolds numbers is clearly seen in the trend shown by the simulations of Bryden et al. (1999). The gap width increases more slowly than predicted by equation(10) and is off by a factor of $`2`$ or more at low viscosity. Again this is consistent with the fact that equation(10) does not appropriately estimate gap widths when the edge of the disk is near the lower order Lindblad resonances. It is difficult to draw firm conclusions from this data since Bryden et al. (1999) do not explore the relation across a uniform (in log) sampling of viscosities. Figure 1 also shows data from our simulations with a more complete sampling of viscosity and these confirm the weak dependence of gap width on Reynolds number.

Thus it appears that while previous approaches can predict when a gap forms they do not capture the full physics of the gaps. In particular some aspect of the scenario in which the torque from waves driven at the OLR is balanced by viscous accretion (as discussed by Bryden et al. 1999) must be incomplete. This is not a trivial point as observational diagnostics from platforms such as SIRTF and ALMA may determine gap widths to less than 10% (Wolf et al., 2002). Since many of the observable systems may be appropriate for the low viscosity limit, understanding the proper scaling there will be critical to using observations of gap widths and disk edges to infer protoplanet properties.

## 3 Numerical Code Description

To investigate the evolution of gaps opened by a protoplanet, we ran a set of $`2`$D hydrodynamical simulations using the code developed by Masset (2002); Masset & Papaloizou (2003). This code is an Eulerian polar grid code with a staggered mesh and an artificial second order viscous pressure to stabilize the shocks (see also Stone & Norman 1992). The hydrocode allows tidal interaction between one or more planets and a 2D non-self-gravitating gaseous disk, and is endowed with a fast advection algorithm that removes the average azimuthal velocity for the Courant timestep limit (Masset, 2000). The simulations are performed in the non-inertial non-rotating frame centered on the primary star (similar to a heliocentric frame). The outer boundary does not allow either inflow or outflow, so it must be located sufficiently far from the planet to ensure that spiral density waves are damped before they reach it. This is facilitated by adopting a logarithmic grid in radius. The grid inner boundary only allows material to escape so that the disk material may be accreted on to the primary star.

The code is scaled so that the unit length is $`1`$ AU which is the initial position of the protoplanet, the unit mass is $`1M_{}`$, which is the mass of the central star, and time is given in units of one orbital period (one year) for an object in a circular orbit at the initial location of the planet.

The runs we are presenting here are made with a resolution of $`N_r=150`$ and $`N_\theta =450`$. The grid physical size spans radii between $`0.25`$ AU and $`6.0`$AU. The viscosity is parametrized with the Reynolds number $`r^2\mathrm{\Omega }/\nu `$, where $``$ is assumed to be constant with radius. The disk aspect ratio $`H/r=0.04`$ is uniform and constant where $`H`$ is the vertical scale height of the disk. The sound speed of the gas is set from the disk aspect ratio. The planet was initially embedded in a disk with density profile $`\mathrm{\Sigma }(r)=10^5\left(\frac{r}{R_{min}}\right)^{1.5}`$ where $`R_{min}`$ is the radius of the inner edge of the grid. Table 1 shows the characteristics of each run we will discuss here.

## 4 Results

Figure 2 show snapshots of the density for a typical simulation, namely $`q=10^3`$ (corresponding to a Jupiter mass planet) and a Reynolds number $`=2\times 10^6`$. From top to bottom and left to right we show the density distribution at times $`t=20,50,100,200,500,1000,2000`$ and 3000 (in units of the planet’s orbital period). We see that a gap is cleared relatively quickly during the first hundred rotation periods. During this time period a prominent 2-armed structure is seen in the density distribution. At later times, the density near the planet drops and the amplitude of the spiral density waves decreases. An asymptotic state is approached on longer timescales when the edge of the disk is between the $`m=2`$ and $`m=1`$ OLRs. We see that the gas interior to the planet slowly accretes onto the star while that exterior to the planet slowly accretes inward piling up near the edge of the disk, outside the planet.

The shape of the disk edge at latter times is better seen on Figure 3a which shows the density profile of the inner region of the disk for the same sampled times. Figure 3a also shows 3 radial functions that can be used as metrics to measure $`\mathrm{\Delta }`$, the “outer” width of the gap or the distance between the planet and the gap outer edge.

### 4.1 Measuring the distance between the planet and the edge of the disk

To measure the distance between the planet and the edge of the disk from a numerical simulation we must first consider procedures for defining a point of measurement for the disk edge. As seen in Figure 3a, the gaps are not well approximated by a square well. Figure 3a illustrates three procedures for defining a location to measure the disk edge from a radial density profile. In each case we choose a stationary radial function and define $`\mathrm{\Delta }`$ based on the location where the radial density profile crosses this function. Our first case, defined by the upper line in Figure 3, is the initial unperturbed density profile. For the second case we use one half of the initial density profile to define our points of measurement. Our third case (represented by the horizontal line in Figure 3) is defined with an absolute density of $`\mathrm{\Sigma }=10^6`$ (independent of radius) which is approximately an order of magnitude below the initial disk density at $`r=r_p`$.

We now compare the evolution of $`\mathrm{\Delta }`$ based on these measurements. From Figure 3b we see that the second and third definitions evolve similarly, but the first one, which refers to the initial density profile, continues to increases on long timescales to a greater extent than the other measurements. Because accretion continues from the outer parts of the disk, the gas density at the outer boundary of the gap continues to increase on long timescales. It is likely that this slow accretion governs the longterm behavior at the “top” of the density profile. To minimize the sensitivity to the exact shape of the edge profile we subsequently choose to measure the distance between the planet and the edge of the disk using our third measurement case (referring to an absolute density of $`\mathrm{\Sigma }=10^6`$). Because there are no significant differences between the second and third measurements, we have chosen to subsequently measure the simpler of the two.

### 4.2 Evolution of the disk edge

Figure 4 shows the evolution of $`\mathrm{\Delta }`$, the distance between the planet and the disk edge, and the density at the location of the planet for three simulations (runs #6, #8 and #11 listed Table 1) which share either the same planet mass ratio or viscosity. The starred points show the depth and width for a planet mass ratio of $`q=2\times 10^3`$ in a disk with $`=2\times 10^6`$. We refer to this run as the “reference run”. The crosses show a run that has the same Reynolds number but a planet mass ratio $`q=10^3`$, half that of the reference run. The last simulation (shown with diamond points) has the same planet mass ratio as the the reference run but a higher disk viscosity ($`=2\times 10^5`$).

Our initial density distribution, that of a smooth power law, is not an equilibrium state in the presence of a planet. As previous numerical studies have illustrated, the planet accretes gas and tidally induces spiral density waves that repel the gas from the planet. Figure 4 shows that two regimes govern the evolution of the depth and width of the gap. On short timescales ($`t<200`$ orbital periods) the clearing of gas near the planet is rapid. The depth of the gap is almost independent of the planet mass and viscosity, however the distance between the planet and the gap outer boundary of the disk is dependent upon both. After $`t200`$ orbital periods, both the density and gap width evolve much more slowly, as if approaching a steady state.

A careful inspection of Figure 4 shows that the initial evolution of the simulations are very similar, and nearly independent of the disk viscosity. This suggests that at early times the angular momentum flux is dominated by spiral density waves driven by the planet, and that the viscous inflow rate is relatively minor. The clearing of the initial gap should take place on a timescale that is set by how fast the spiral density waves can push mass away from the planet. The torque on the disk can be estimated from the sum over each resonance outside the edge of the disk. If the density drops near the planet, then the Lindblad resonances nearest the planet will not significantly contribute to the torque on the disk. In other words, because equation (3) depends on the density at the resonance, only resonances outside the gap in the disk will exert a significant torque on the disk.

At the beginning of the simulation the gas disk extends all the way to the Roche lobe of the planet and the clearing timescale is set by a summation over high order (large $`m`$) Lindblad resonances. However as gas is cleared away from the planet, the density drops at the highest order resonances and the lower order Lindblad resonances take over. Since these aren’t as strong, this clearing should go at a slower rate. So we expect that much of the clearing takes place under the influence of only the last few resonances. From equation(7) we expect evolution on a timescale in planetary rotation periods

$$t_{clear}\frac{r_p^4\mathrm{\Omega }_p^2\mathrm{\Sigma }}{2\pi _mT_m}250\left(\frac{q}{10^3}\right)^2$$
(12)
where we have added the contributions from the $`m=`$2,3 and 4 Lindblad resonances.

When the gap opens sufficiently that the higher order resonances are evacuated and the $`m=2`$ OLR resonance dominates, the timescale could be an order of magnitude longer. Nevertheless, the above timescale is approximately what we see for the first phase of evolution from the simulations presented in Figure 4, and is consistent with this first phase of evolution being independent of the disk viscosity. Since the timescale for clearing the gas away from the planet depends primarily on the planet mass, we expect that the simulations with massive planets will initially clear material faster than those with lower mass planets. This is primarily seen in the evolution of the disk edge (Figure 4b) and not in the depth of the gap (Figure 4a), though the depth of the gap in the lower planet mass simulation (crosses) is higher than that of the higher planet mass simulation at the same viscosity (stars).

On longer timescales viscous behavior becomes important. For the higher viscosity simulations we expect the viscosity becomes important earlier on, consistent with the quicker transition toward the regime of slower evolution seen in the higher viscosity simulation (diamonds) in both density and width plots of Figure 4. Except for the lowest viscosity simulations considered in this paper, the transition to the longer slower evolution takes place near $`t=200`$ orbital periods. So measurements taken after this time can be considered to represent quantities measured when the disk is slowly evolving. Because accretion continues from the outer disk, gas continues to pile up at the edge of the disk. Consequently these simulations never reach a steady state, but on long timescales slowly evolve on the viscous timescale of the outer disk.

The long term accretion complicates the development of an understanding of the scaling behavior of the edge profile on planet mass and disk viscosity. To enable us to explicate the dependence of the disk edge profile on planet mass and disk viscosity we will subsequently compare gap properties measured at different times and discuss the scaling of the properties of disk edges keeping in mind that our simulations show evolution on two timescales.

### 4.3 Influence of the disk’s viscosity

Figure 5 shows the depth and distance to the edge of the disk for a simulation with planet mass ratio $`q=2\times 10^3`$ as function of the Reynolds number. The upper curve (diamond points) is at $`t=200`$, i.e. in the first phase of evolution. The lower curve (starred points) is at $`t=1000`$, and shows the second phase of evolution. At latter times when the viscous timescale becomes important, we see that the depth is more strongly influenced by the disk viscosity than during the first phase. The depth at $`t=200`$ is independent of viscosity for $`10^6`$ because the secondary timescale evolution has not yet completely taken over (see Figure 4).

Although the gas density near the planet is strongly dependent on the viscosity at later times, the distance between the planet and the edge of the disk is only weakly dependent upon the disk viscosity (see Figure 5b). A variation of 2 orders of magnitude in the disk viscosity causes only a 30% change in the distance between the planet and disk edge, but a similar scale (2 orders of magnitude) variation in the gas density. The strong dependence of the gas density near the planet on viscosity, and weak dependence of the actual location of the disk edge on the disk viscosity will be discussed in section 6 when we consider modifications to the previous theory.

### 4.4 Influence of the protoplanet’s mass

Figure 6 shows the behavior of the gap and disk edge created in a disk of viscosity corresponding to $`=2\times 10^6`$ as function of the planet mass ratio. We plot results for simulations carried out with planet mass ratios between $`q=10^4`$ and $`2\times 10^3`$. The two lines shown in both figure panels are similar so we can study the dependence of the gap density and width on the planet mass ratio without being sensitive to the time at which the measurements were taken.

We see from Figure 6 that both the depth and the distance between the planet and disk edge vary significantly as function of the planet mass even though we have only run simulations over a relatively small range in planet mass ratio (particularly when we compare to the larger range of Reynolds numbers that we explored). Since a gap will not open at the lowest planet masses and higher ones exceed the masses of most known extrasolar planets, we cannot significantly extend the range of planet masses considered. In the regime plotted in Figure 6 we see that a variation in the planet mass ratio of a factor of 4 causes approximately an order of magnitude a change in the gas density near the planet, and a 40-50% change in the distance between the planet and the disk edge. As we previously found for the dependence on disk viscosity, the distance between the planet and the disk edge is less strongly dependent upon the planet mass ratio than is the gas density near the planet. Both the depth and distance between the planet and disk edge are more strongly dependent upon the planet mass ratio than on the disk Reynolds number. In the following section we will attempt to account for the trends we observe in the long timescale evolution of these simulations.

## 5 Comparisons to Theory

As suggested by Bryden et al. (1999), when mass has been removed from the region nearest the planet and a steep density gradient is present, diffusion or viscous spreading may become important and determine the shape of the edge of the disk (see Figure 7). The basic equation that describes the evolution of a Keplerian disk due to viscous processes is

$$\frac{\mathrm{\Sigma }}{t}=\frac{3}{r}\frac{}{r}\left[r^{1/2}\frac{}{r}\left(\nu \mathrm{\Sigma }r^{1/2}\right)\right]$$
(13)
(Pringle, 1981). This equation arises by considering the Navier-Stokes equation and conservation of mass. The azimuthal component of the Navier-Stokes equation, which addresses angular momentum transport, is used to estimate the radial velocity component. This, when inserted into the equation for conservation of mass, yields the above equation. When the density $`\mathrm{\Sigma }`$ is nearly flat, the angular momentum flux or torque through a radius $`r`$ is given by equation(8). This is what is usually meant by “accretion,” as we shall see. When the radial gradient of $`\mathrm{\Sigma }`$ is very large however, there is an additional diffusive term that drives mass flow.

The previous equation can be modified to include mass and angular momentum flow due to dissipation of spiral density waves driven by the planet by adding a term,

$$\frac{\mathrm{\Sigma }}{t}=\frac{3}{r}\frac{}{r}\left[r^{1/2}\frac{}{r}\left(\nu \mathrm{\Sigma }r^{1/2}\right)+\frac{1}{3\pi r\mathrm{\Omega }}\underset{m}{}\frac{F_m}{r}\right]$$
(14)
where $`\frac{F_m}{r}`$ corresponds to the torque density exerted in the disk due to dissipation of spiral density waves driven by the planet at the $`m`$th Lindblad resonance, and there is a sum of terms due to the contributing resonances. Here we have followed the derivation by Takeuchi et al. (1996).

Once a prominent gap is opened we must distinguish between the physics occurring within the gap, (particularly along its steep outer edge), and that occurring outside the gap edge at larger radii. Outside of the gap the density gradient is shallow, the diffusion rate will be small and we expect that the mass inflow due to accretion will dominate and will be balanced by mass outflow resulting from the spiral density waves. However, in the steep edge itself, diffusion should be important. In this region it is mass inflow due to diffusion that should be balanced by the mass outflow caused by dissipation of the spiral density waves. Thus we can divide the problem into two regimes: accretion dominated and diffusion dominated (see Figure 7). This division can be directly seen in Figure 3a showing the density profile of a single simulation. Past the edge of the disk ($`r1.3`$) the density gradient is small, but inside this radius, the density drops by orders of magnitude and the density gradient is large.

Unless the planet mass is sufficiently low and viscosity sufficiently high that a gap is only barely opened, we expect that edge of the disk will be steep. The resonance most distant from the planet at highest density $`\mathrm{\Sigma }`$ should dominate the shape of the outer part of the edge because it will be most effective at driving spiral density waves. We refer to quantities such as radius and density past the edge of the disk as $`r_e`$ and $`\mathrm{\Sigma }_e`$ and those at the $`m`$-th OLR resonance as $`r_m`$ and $`\mathrm{\Sigma }_m`$ (see Figure 7).

In the limit that the density is not varying quickly with time, ($`\frac{\mathrm{\Sigma }}{t}0`$), the previous equation can be written

$$3\nu \mathrm{\Sigma }+3\nu r\frac{\mathrm{\Sigma }}{r}\frac{1}{3\pi r\mathrm{\Omega }}\underset{m}{}\frac{F_m}{r}+\mathrm{constant}$$
(15)
where the first term on the left is what is typically known as accretion<sup>1</sup><sup>1</sup>1 The first term is $`3\nu r\mathrm{\Sigma }/2`$ when $`\nu `$ is independent of $`r`$ leading to equation 8, however here we have assumed that the Reynolds number is constant, consistent with an ”$`\alpha `$-disk”. and the second can be described as diffusion. Here $`F_m`$ is the angular momentum flux transfered by the waves and $`\frac{F_m}{r}`$ is the torque density exerted on the disk due to the dissipation of the waves in the disk (Takeuchi et al., 1996). Using the definition for the Reynolds number, we can rewrite the previous equation as

$$\frac{(r\mathrm{\Sigma })}{r}\frac{}{3\pi GM_{}}\underset{m}{}\frac{F_m}{R}.$$
(16)
Integrating the previous equation we find

$$r\mathrm{\Sigma }\underset{m}{}\frac{F_m}{3\pi GM_{}}+\mathrm{constant}$$
(17)
which implies that the shape of the disk edge is entirely determined by the way that the waves are dissipated. The steep edges we see in the simulations imply that the dissipation occurs very rapidly, as previously pointed out by Takeuchi et al. (1996).

The angular momentum flux transferred by the spiral density waves $`F_m`$ is equal to that given by equation (3) at the location of the resonance $`r_m`$ and then decays because of dissipation. We assume that the spiral density waves are dissipated over a length scale $`r_d`$ (referred to as a damping length by Takeuchi et al. 1996), so $`F_mT_m(1e^{(rr_m)/r_d})`$ for $`r>r_m`$ and $`T_m`$ is given by equation (7) with quantities defined at the resonance. We associate the edge of the disk as that point where the dissipation of the density waves ceases to be important; $`r_er_m+r_d`$ where $`r_m`$, in this case, is the outermost resonance that still lies within the disk edge.

Considering only one resonance we can rewrite the previous equation as

$$r\mathrm{\Sigma }r_m\mathrm{\Sigma }_m+\frac{\mathrm{\Sigma }_mf_mr_p}{3\pi }\left(\frac{q}{10^3}\right)^2\left(1e^{(rr_m)/r_d}\right)$$
(18)
for $`r>r_m`$. We expand this to include all resonances.

$$r\mathrm{\Sigma }r_p\mathrm{\Sigma }_p+\underset{m}{}\frac{\mathrm{\Sigma }_mf_mr_p}{3\pi }\left(\frac{q}{10^3}\right)^2a_m(r)$$
(19)
where

$`a_m(r)=`$ $`0`$ $`\mathrm{for}r<r_m`$ (20)
$`1e^{(rr_m)/r_d}`$ $`\mathrm{for}r>r_m.`$
First let us consider the outermost resonance that still lies within the edge of the disk. Since the rightmost term of this equation is large for $`r>r_m+r_d`$, we can estimate the density ratio at the edge compared to that at the resonance

$$\frac{\mathrm{\Sigma }_e}{\mathrm{\Sigma }_m}\frac{f_m}{3\pi }\left(\frac{q}{10^3}\right)^2.$$
(21)
We find that the drop in density near the disk edge is dependent on the square of the planet mass ratio and the Reynolds number. This is consistent with our numerical exploration which found that the gap depths were strongly dependent on both quantities and more strongly dependent on the planet mass ratio than the disk viscosity (section 4.3, 4.4). However, the distance from the outermost resonance to the peak density, $`\mathrm{\Sigma }_e`$, only depends on the scale over which the dissipation takes place, $`r_d`$. Using the WKB approximation, Takeuchi et al. (1996) found that the dissipation scale length was primarily dependent on the Mach number of the disk, but also dependent on the disk viscosity. The weak dependence in our simulations we found for the distance of planet to the disk edge implies that the dissipation length is not strongly dependent on the viscosity. This may in part be due to the formation of shocks (e.g., Ward 1997).

Equation (21) allows us to estimate the density at the radius of the dominant Lindblad resonance from the following parameters; the resonance parameter $`f_m`$, the planet mass ratio, $`q`$, and the Reynolds number, $``$. We consider the simulation displayed in Figures (5, 6) with planet mass ratio $`q=2\times 10^3`$ and $`=2\times 10^6`$. Because the gap outer edge is nearest the $`m=2`$ OLR, we expect that this resonance, with $`f_2=5\times 10^5`$, balances the viscous inflow. From these parameters and using equation (21) we estimate $`\frac{\mathrm{\Sigma }_e}{\mathrm{\Sigma }_2}42`$. This ratio and the scale length $`r_d`$ determines the slope in the outermost brightest and most easily detectable region of the disk edge. While the ratio $`\frac{\mathrm{\Sigma }_e}{\mathrm{\Sigma }_2}`$ is greater than 1, it is not high enough density ratio to explain the orders of magnitude lower densities seen the planet that we measured in the simulations (see Figures 5,6). Consequently we conjecture that more than one resonance must affect the the disk edge profile.

To consider the role of the different resonances we look more closely at the density profile shown in figure 3a. Equation (21) allows us to estimate the change in density between the outermost resonances. For $`q=10^3`$, $`=2\times 10^6`$, we expect $`\frac{\mathrm{\Sigma }_e}{\mathrm{\Sigma }_2}11`$, accounting for order of magnitude drop in density between $`r=1.3`$ (the location of the $`m=2`$ OLR) and $`r1.5`$ (the location of peak density) as is consistant with what is seen in figure 3a. Between the $`m=3`$ OLR (at $`r=1.2`$) and the $`m=2`$ OLR (at $`r=1.3`$) we expect the density to change by a factor of $`\frac{f_3}{3\pi }\left(\frac{q}{10^3}\right)^2=42`$ which again is consistant with for the drop in density between $`r=1.3`$ and $`r=1.2`$ of another order of magnitude.

Thus it appears we can divide the disk edge into two regimes, one it which the distance between the resonances is larger than the dissipation scale length $`r_d`$ and another where the distance between the resonances is smaller than $`r_d`$. Note that the distance between the resonances increases as a function of distance from the planet and is inversely dependent on $`m`$. When the distance between the resonances is larger than the dissipation scale length then the waves can fully dissipate their angular momentum and the density increases by a factor of $`\frac{f_m}{3\pi }\left(\frac{q}{10^3}\right)^2`$. Each resonance causes an increase in density by a multiplicative factor which depends on the Reynolds number, the square of the planet mass and $`f_m`$.

We now consider the regime nearer the planet where the distance between the resonances is smaller than the dissipation scale length. In this regime we the waves do fully dissipate their angular momentum before the next resonace is encountered. We conjecture that the effect of the resonances is additive in this case rather than multiplicative. Note here that when the resonance Mach number $`\frac{r\mathrm{\Omega }}{mc}<1`$ the resonances can no longer efficiently drive spiral density waves. Thus there is a cutoff in $`m`$ and a lower limit on the density near the planet.

Because the higher order resonances only additively affect the density profile, the ratio of the density near the planet compared to that outside the disk edge is primarily determined by the number of resonances in the outer regime, where the distance between the resonances is smaller than the dissipation scale length and the density increases by a factor of $`\frac{f_m}{3\pi }\left(\frac{q}{10^3}\right)^2`$ between each resonance. The high power dependence of the ratio of the density between the planet and disk edge on $`q`$ and on $``$ (that we remarked on in sections 4.3 and 4.4; see Figures 5 and 6) can be explained if more than one resonance is in this regime, and if $`r_d`$ is not strongly dependent on the planet mass or disk Reynolds number.

In contrast, we expect the location of the disk edge to be set by the outermost resonance capable of driving significant density waves and the dissipation scale length $`r_d`$. From equation (19) we expect that the outermost resonance is the lowest $`m`$ for which $`\frac{f_m}{3\pi }\left(\frac{q}{10^3}\right)^21`$. This limit is approximately consistent with the distance between the planet and disk edge dropping below 0.3 for $`q10^3`$ seen in figure 6b. In section 4.3 and 4.4 we found that the distance between the planet and disk edge was only weakly dependent on the planet mass and disk viscosity. Once the outmost resonance is set, equation(19) suggests that that this distance only depends on the dissipation scale length. If this is not strongly dependent on the disk visocity then we can account for weak dependence of this distance on planet mass and Reynolds number.

## 6 Summary and Discussion

In this paper we have presented a series of 2D numerical hydrodynamical simulations of planets embedded in Keplerian disks. We have covered a different area of parameter space compared to previous works, and have concentrated on lower viscosities or higher Reynolds number disks with the expectation that this regime would represent a more extended time after planets have formed and have ceased to accrete gas from the protostellar disk. As previous simulations have shown, the planet can open a gap in the disk. On longer timescales we find the width of the gap depends only weakly on the mass of the planet and disk viscosity. The depth of the gap and slope of the outside edge of the gap is, however, comparatively more strongly dependent on both quantities, and a stronger function of planet mass ratio than of viscosity. When $`q^20.3`$, the edge of the disk can be found between the $`m=2`$ OLR and the $`m=1`$ OLR. This is unexpected since the planet can only resonantly drive spiral density waves at the locations of resonances.

We account for the phenomena seen in the simulations with a scenario model which balances the viscous transport of angular momentum against that arising from the dissipation of spiral density waves induced tidally by the planet. The disk approaches a steady state where outside the gap the inward mass flux from viscous accretion is balanced by outward flux due to the dissipation of planet driven spiral density waves. Within the steep gap edge, however, it is diffusion that balances the outward angular momentum flux from the dissipation of the spiral density waves. These waves are excited inefficiently because of the extremely low gas density near the planet. The profile shape is determined by the scale length over which the spiral density waves are dissipated, $`r_d`$.

From our simplified model we conjecture that the edge can be divided into two regimes: where the distance between resonances is greater than $`r_d`$ and where this distance is smaller than the dissipation scale length. In the outermost regime each resonance multiplicatively increases the density as a function of distance from the planet. The multiplicative factor depends on the product of the Reynolds number and square of the planet mass. This prediction matches the stronger dependence of gap properties on planet mass ratio then on Reynolds number that we also observed in the simulations. Finally the distance between the planet and the disk edge is set by the outermost resonance capable of driving strong density waves; the $`m`$-th OLR for which $`\frac{f_m}{3\pi }\left(\frac{q}{10^3}\right)^21`$ and by the dissipation scale length. The weak dependence we saw for this distance on the planet mass and Reynolds number suggests that the dissipation scale length is not strongly dependent on the disk viscosity. Once the outermost resonance is determined, the slope of outer, most easily detected part of the disk edge depends on $`r_d`$ and the ratio $`\frac{f_m}{3\pi }\left(\frac{q}{10^3}\right)^2`$.

Rice et al. (2003) noted that the slope of a disk edge maintained by a planet is dependent upon the planet mass and suggested that this dependence could be used to constrain the type of planet that could be present from fits to the spectral energy distribution (as from GM Aurigae). The scaling we have developed confirms this work and furthermore suggests that the slopes of disks edges maintained by planets are approximately proportional to the product of the planet mass ratio squared and the Reynolds number.

The large drop in density in the disk edges predicted by our simulations implies that the opacity of the gas near the planet is much reduced. Because disk edges are illuminated by the central star, and could be optically thin, the spectral energy distribution predicted from the simulations may be sensitive to the slope of the edge. As proposed by Rice et al. (2003), its possible that circumstellar disk edges may be used to constrain the masses of planets which maintain them. We note that because the outer disk slope depends on both $`r_d`$ and a density ratio which is set by the outermost resonance, there could be more than planet mass setting observed outer edge slopes. Future work could focus on how the scale length $`r_d`$ affects the slope and so might discriminate between possible degenerate solutions. One interesting possibility is that the changes in gas density slope can be detected. The ratio of the slopes on either side of the transition point is set by the resonance at the transition region and the distance from the transition point to the density peak is set by the dissipation scale length. If such measurements became possible then the degeneracies could be lifted and the planet mass tightly constrained from the disk edge density profile.

Our simulation results clearly point to the need for a deeper understanding of the physics of gaps and disk edges. Simplistic models such as explored in section 2, that neglect diffusion and the density gradient in the torque balance estimate, would predict a stronger scaling on the gap size with $`q`$ and $`R`$ than we have measured from our simulations. The diffusive model we discuss here is preliminary because it is used to approximate long timescale behavior seen in the simulations which never actually achieve a steady state. We have neglected variations in the angular rotation rate caused by the pressure differential in the steep edge of the disk. We estimate however that the location of the resonances are not significantly affected by this pressure differential. Our model also assumes that spiral density waves are driven at the location of the resonance and we have not considered models that drive density waves at significant distances from Lindblad resonances (e.g., as exhibited in the simulations of Savonije et al. 1994). The ability of a planet to drive density waves off resonance depends on the sound speed of the disk (and possibly on the boundary conditions) and it is more likely this would occur in low Mach number disks. In this situation we expect that the gas density near the resonances would be reduced to even a larger extent to balance viscous inflow in the high planet mass, low viscosity simulations. The importance of this driving mechanism and the resulting affect on the disk edge can be explored by carrying out simulations of thicker disks. By affecting the radial range over which the spiral density waves are launched, the thickness of the disk may influence the shape of the edge profile. The simulations of Kley (1999) suggest that thicker disks have shallower disk edges. This sensitivity, which we have not yet explored, would increase the number of parameters required to describe a disk edge. Here we have adopted an exponential form for the dissipation of the waves with radius. However the sensitivity of the dissipation scale length to sound speed and viscosity has not yet been explored beyond that discussed by Takeuchi et al. (1996) and the sensitivity of the disk edge profile to these parameters could be more carefully explored. Our work has been restricted to the case of a planet in a circular orbit, however the interaction between a planet and a disk may induce eccentricity in the planet (Goldreich & Sari, 2003). In this case, because of the increased number of resonances at which spiral density waves can be driven, it should be even more difficult to predict the shape of the disk edge profile.

Because mass continues to pile up outside the outer edge of the disk near a planet, a steady state will not be reached, increasing the difficult of predicting the edge profile. The density near the planet will increase until the planet begins to migrate. As disks age and the viscosity drops, the timescale for planetary migration times will also increase, so observed protoplanetary and debris disks may be slowly evolving rather than steady state. Our partial understanding based on the long timescale behavior of these systems gained from this study may be relevant toward the interpretation of disk edges observed in forthcoming observations. Since the edges of disks receive and process radiation differently than disk surfaces, observations may be particularly sensitive to the 3 dimensional shapes of these edges. Future work should also concentrate on developing 3D models for disk edges in the physical regime studied here.

Finally, in a model where planetary formation is sequential, it is tempting to speculate that the gas that piles up on long timescales outside the $`m=2`$ OLR may induce formation of an additional planet (Lin & Papaloizou 1993). Our study suggests that the secondary planet would seldom form locked in the 3:2 mean motion resonance with the primary as is Pluto with Neptune (at $`r1.3`$, the 3:2 mean motion resonance is equivalent to the $`m=2`$ OLR) but would be more likely to form locked in a more distant resonance such as the 2:1 mean-motion resonance (at $`r=1.6`$, as are a number of the extrasolar planets) or in between, as is Saturn with Jupiter.

We thank F. Masset for his help and for providing the code. We also thank an unamed referee for considerable help in making the theory section of the paper stronger by pointing out a number key conceptual issues. Support for this work was provided by NSF grants AST-9702484, AST-0098442, NASA grant NAG5-8428, DOE grant DE-FG02-00ER54600, and the Laboratory for Laser Energetics. This material is based upon work supported by the National Aeronautics and Space Admininstration under Grant No. NAGxxxx issued through the Origins of Solar Systems Program. This research was supported in part by the National Science Foundation to the KITP under Grant No. PHY99-07949.


text:
# Chapter 1 A Study of the Ground Level Event of April 15, 2001 with GRAND

## Chapter 1 A Study of the Ground Level Event of April 15, 2001 with GRAND

### Abstract

Project GRAND, a proportional wire chamber array, is used to examine the increased counting rate of single muons during the solar Ground Level Event (GLE) of April 15, 2001. Data are presented and compared to our background rates and to that of neutron monitor stations. This GLE is seen with a statistical significance of 6.1$`\sigma `$.

### 1. Introduction

On April 15, 2001 the IPS Radio and Space Services website reported a rare X14 ($`14\times 10^4W/m^2`$ at one AU from the sun) solar flare with a maximum X-ray intensity at 13.83 UT. The flare originated on the sun at coordinates S20 W85. During the period from 1976 - 1998, only sixteen flares of magnitude X10 or greater were reported; observations on some of these high energy flares are given in .

During the peak time of the flare, the sun was at an altitude of $`35^{}`$ and an azimuth of $`110^{}`$ in a local coordinate system. Even though the sun is near the minimum detectable elevation angle (27), the interplanetary magnetic field (IMF) and the magnetic field of the earth alter the direction of charged particles from the sun allowing their charged secondaries to be more easily detected by GRAND. Muon telescopes are typically sensitive to higher primary energies than neutron monitors.

Charged particles emitted from the sun spiral around the magnetic field lines which themselves curve in an Archimedes spiral. Therefore the GLE associated with a particular X-ray flare should occur at a slightly later time than the X-ray event (due to the longer path length and slower velocity) and at a different angle than the sun, thus providing better detector acceptance even though a direct line to the sun is at a low elevation angle (see also ).

Further details on the observation of this GLE can be found in .

### 2. Data and Observations

Project GRAND is an array of 64 proportional wire stations located at $`41.7^{}`$ N and $`86.2^{}`$ W at an altitude of 220 m above sea level. The array is used to detect ground level charged particles produced from primaries undergoing hadronic interactions. Each station contains eight 1.29 m2 proportional wire chamber (PWC) planes. The position sensitive PWC planes allow the direction of a muon track to be measured to within 0.26, on average, in each of two projected planes. A 50 mm thick steel plate is located above the bottom two PWC planes. Muon tracks can thus be distinguished from electrons or hadrons. The array collects data at a rate of $``$2400 muons per second. Additional details are on the webpage: http://www.nd.edu/$``$grand.

GRAND recorded a continuous data file on magnetic tape containing information from April 14, 2001 at 8:40 UT through April 16 at 21:30 UT, providing background information before and after the flare time. The counting rate from 14.5 hours to 24.5 hours UT was examined in 0.1 hour bins. Deviations (r.m.s.) from average counting rate were measured individually for each station of the array. In order to prevent spurious individual station variations from giving a possible signal in the sum-of-all-stations, one-third of the stations with the highest individual r.m.s. deviations were eliminated from consideration leaving 39 stations (an extremely severe cut to ensure with great confidence that erratic station behavior could not cause a false signal).

The six bins of time from 14.0 to 14.6 hours UT were used as a signal time; the bins from 9.5 to 19.5 hours UT (minus the signal time) were used as background. Since it is known that the muon rate depends slightly on the time of day  (due to small solar effects in air pressure, temperature, and in the IMF), the counting rate for the background time period was fit to a quadratic curve. The muon counting rates for April 14 and April 16 (the preceeding and the next day from the signal) were examined and it was verified that the backgrounds from these two days are similar to the curve for April 15; the average of four years of data also shows a similar shape. The data as percentage above background are shown in Figure 1 with error bars. For comparison, the data from neutron monitors at Newark  and Oulu  during this GLE are also shown. Our muon data rate has not been corrected for changes in pressure. However, the pressure was recorded every half hour and found to be constant for the entire interval of time analyzed to within $`\pm `$0.15% and thus could not cause the excess counting rate observed. The amount of signal above background yields a GLE signal with a significance of 6.1$`\sigma `$ (where $`\sigma `$ is the total error on the signal and background).

Data from the Newark , and Oulu  Neutron Monitor Stations are compared in Figure 1, the data for Climax  are very close to that of Newark and for clarity are not shown. The energy thresholds for the Climax, Newark, and Oulu neutron monitor stations are influenced by their vertical geomagnetic cutoff rigidity (3.0, 2.1, and 0.8 GV respectively) and overburden of air (their corresponding heights above sea level are 3400, 50, and 15 m). The corresponding numbers for GRAND are a geomagnetic cutoff rigidity of 1.9 GV and a height of 220 m above sea level. However, our primary energy is also influenced by the different mechanisms which produce muons (rather than neutrons) which subsequently reach detection level. While looking at single muon tracks at ground level, our threshold primary proton energy is 3 - 10 GeV , depending on the primary spectral index. Thus our peak is smaller in amplitude due to lower fluxes higher primary energies. The Newark signal has a half-height time at 14.21 hours, the Climax signal at 14.16 hours, and Oulu at 14.27 hours (where the half-height time is the time at which the data has risen above background to half of the highest measured data point of the GLE). A similar calculation for our half-height time is 14.09 hours which is $``$ 0.1 hour earlier (the earlier time might be expected on the basis of higher primary energies compared to the neutron monitors  with faster protons traveling at smaller pitch angles yielding shorter path lengths. The width (full width at half-of-maximum above background) of the Climax, Newark, and Oulu peaks are 1.27, 1.24, and 1.61 hours in contrast to our value of 0.48 hour. Our peak is thus $``$0.9 hours shorter in duration (typically, higher energies have a shorter duration than lower energies as is seen in satellite data for protons at lower energies.

### 3. Conclusions

Project GRAND sees a ground level signal with a significance of 6.1$`\sigma `$ when examining the secondary muon counting rate at ground level between 14.0 and 14.6 hours UT on April 15, 2001. This signal is obtained with no restrictions on the angle of the muons, consistent with protons originating at the surface of the sun and accelerated during the time of the X14 X-ray flare. Our data at $`\pm 5^{}`$ angles from the direction of the sun show no increase; therefore our detection is not due to gamma rays or neutrons from the sun. Our GLE detection occurs $``$0.1 hour earlier than Climax, Newark, or Oulu neutron monitor signals and is shorter in duration which might be expected considering the higher energies involved.

The authors wish to thank M. Dunford, J. Vermedahl, and M. López del Puerto. Thanks to L. Miroshnichenko for preliminary intensity calculations and M. Shea and D. Smart for calculations of our asymptotic direction as a function of rigidity; The GOES science team, J. Clem and R. Pyle (Bartol Neutron Monitor Program funded through NSF Grant ATM-0000315), C. Lopate (Climax Neutron Monitor operated by The University of Chicago and funded through NSF grant ATM-9912341), and the Oulu Neutron Monitor for the use of their data. GRAND is funded through the University of Notre Dame and private grants.

### 4. References

1. Bieber J. et al. 2002, ApJ 567, 622

2. Clem J. and Dorman L. 2000, Space Science Review 93, 335

3. Clem J. and Pyle R. 2002, Private Communication

4. Falcone A., Ryan J., 1999, Astroparticle Physics 11, 283

5. Lopate C. 2002, Private Communication

6. Miroshnichenko L. 2001, Solar Cosmic Rays, (Kluwer, Norwell, MA)

7. Miroshnichenko L. 2002, Private Communication

8. Munakata K. et al. 2001, Proceedings of the 27th International Cosmic Ray Conferece (Hamburg), 3494

9. Oulu Neutron Monitor Website 2002, http://spaceweb.oulu.fi/projects/crs

10. Poirier J., Roesler S., Fasso A. 2002, Astroparticle Physics 17, 441

11. Poirier J., D’Andrea C. 2001, Proceedings of the 27th International Cosmic Ray Conferece (Hamburg), 3923

12. Poirier J., D’Andrea C. 2002, Journal of Geophysical Research 107, 1376

13. Shea M. and Smart D. 2002, Private Communication

14. Swinson D. B., Shea M. A. 1990, Geophysical Research Letters 8, 1073