Initial conditions

• We consider the standard model of protoplanetary disks, where surface density is a power low. We solve the Euler equations, that is to say fluid motion equations with no viscosity. We assume the pressure to suit the perfect gaz low.

Physical model

We use power laws for the physical quantities

• Radius: R = 1.5e13*r cm, so r is a number of AU

• Surface density: Rho = 1700 r^-3/2 g.cm-2, sigref = 1700

• Pressure: P = 16.9107e12 r^-2 Barye.AU, pref = 16.9107e12

• P & Rho*T, with T = 280 r^-1/2 K, and gamma of the gaz is 1.4

• Tangential velocity: Vg = 2.97396e6 r^-1/2 cm.s-1, vref = 2.97396e6, vref² = G.Mstar/(1 AU)

• Radial velocity: Ug = 0 * vref cm.s-1

Numerical disk initialisation

Our numerical code computes the normalized quantities: rosg, pg, ug, vg

• Radial extension: rin = 5, rout = 10
• Linear polar coordinates: r(i+1) = r(i) + dr , theta(j+1) = theta(j) + dtheta, with dr = (rout -rin)/nx and dtheta = 2*pi/ny

Setting the initial bump in variables:

• Surface density: rosg = r^-3/2 * bump(r), so Rho = sigref * rosg

• Pressure: pg = pref/(sigref*vref²) r^-2 * bump(r), so P = sigref*vref² * pg

• Radial velocity: ug = 0, so Ug = vref * ug

• Tangential velocity: vg² = 1/r + r/rosg d/dr pg, so Vg = vref * vg (this is why P=sigref*vref²*pg)

  The bump(r) function is: 1 + A exp[- ((r - 7.5)/drw)²], a gaussian bump. We use A = 0.2 and drw = 0.3 AU

Random perturbation initialization

• In order to trigger the instability, we have added 100 random perturbations on the density and the pressure defined above. The procedure is in Fortran 90:

 call random_seed
 call random_number(t1) ! t1, t2, t3 are arrays of 100 reals, that will receive random numbers in the range [0, 1]
 call random_number(t2)
 call random_number(t3)

 do i1 = 1, 100

    r_rand = 7 + t1(i1) * 1D0   ! adds perturbations in the range [7, 8] AU
    theta_rand = t2(i1) * 2.D0*pi  ! adds perturbations in the range [0, 2 pi] 

    do i = ini, nx
       do j = ini, ny

 !------- Definition of a random gaussian bump -------------
          r_pert = 2.0D0 * (r(i) - r_rand)/drw     ! drw is the same as above (0.3 AU in the example)
          theta_pert = 2.0D0 * (theta(j) - theta_rand)/drw
          bump_pert = 1D0 + 1.D-2 * t3(i1) * dexp(- (r_pert**2 + (r(i)*theta_pert)**2) )

 !------- Perturbating the density and pressure -------------
          rosg(i, j) = rosg(i, j)*bump_pert
          pg(i, j) = pg(i, j)*bump_pert

       enddo
    enddo

 enddo