A Numerical Study of Brown Dwarf Formation via Encounters of Protostellar Disks
Sijing Shen, James Wadsley, Tristen Hayfield, Nicholas Ellens
aa r X i v : . [ a s t r o - ph . S R ] S e p Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 30 October 2018 (MN L A TEX style file v2.2)
A Numerical Study of Brown Dwarf Formation viaEncounters of Protostellar Disks
S. Shen, ⋆ J. Wadsley , T. Hayfield and N.Ellens Department of Physics and Astronomy, McMaster University, Main Street West, Hamilton L8S 4M1, Canada
30 October 2018
ABSTRACT
The formation of brown dwarfs (BDs) due to the fragmentation of early-stageproto-stellar disks undergoing pairwise encounters was investigated. High resolutionallowed the use of realistic, flared initial disk models where both the vertical struc-ture and the local Jeans mass were resolved, preventing numerical enhancement orsuppression of fragmentation. The results show that objects with masses ranging fromgiant planets to low mass stars can form during such encounters from disks that werestable before the encounter. The parameter space of initial spin-orbit orientations andthe azimuthal angles for each disk was explored with a series of 48 simulations. Forthe types of interactions studied, an upper limit on the initial Toomre Q value of ∼ M J , with the ma-jority in the BD regime. These often resided in star-BD multiples and in some casesalso formed hierarchical orbiting systems. Most of them had large angular momentaand highly flattened, disk-like shapes. With the inclusion of the appropriate physicsthese could reasonably be expected evolve into proto-BD disks with jets and subse-quent planet formation around the BD. The objects had radii ranging from 0.1 to10 AU at the time of formation. The disk gas was assumed to be locally isothermal,appropriate for the short cooling times in extended proto-stellar disks but not forcondensed objects. An additional case with explicit cooling that reduced to zero foroptically thick gas was simulated to test the extremes of cooling effectiveness and itwas still possible to form objects in this case. Detailed radiative transfer (not studiedhere) is expected to lengthen the internal evolution timescale for these objects so thatthey spend considerable time as puffed-up, prolate ellipsoids, but not to alter the basicresult that proto-BD disks can form during proto-stellar encounters. Observations suggest that brown dwarfs (BDs, 0.013 M ⊙ . m . . M ⊙ ) are very abundant in both the field and starclusters in our galaxy (Chabrier 2002; Mart´ın et al. 2001;Luhman 2000). The origin of brown dwarfs, however, is stillunder debate. It has been suggested by several authors (e.g.,Whitworth et al. 2007; Bate et al. 2002) that the final stel-lar mass may be closely related to the Jeans mass of theparent gas within the molecular cloud. Given that the typ-ical Jeans mass in a typical molecular cloud is of order onesolar mass, this is reasonable for larger, hydrogen-burningstars. However, because BD masses are at least an orderof magnitude smaller, it is difficult for BDs to form by di-rect collapse to pre-BD cores. Several scenarios have beenproposed to circumvent this problem. One approach is to assume that the initial Jeans massof the core is still a good estimate of the final mass from adirect collapse, but that denser or colder regions exist wherethe Jeans mass may approach the BD mass. For example,Padoan & Nordlund (2004) proposed that local density of incores can be enhanced by turbulent flows (“turbulent frag-mentation”), while Elmegreen (1999) suggested that certaintypes of molecular cooling can greatly decrease core temper-atures. In either case, formation of BDs follows the conven-tional star formation path. Hence, all properties associatedwith proto-stars, such as accretion disks and outflows, areexpected.An alternative picture, first proposed byReipurth & Clarke (2001), suggests that the becausethe Jeans mass continuously decreases during the collapse, c (cid:13) S. Shen et al. fragmentation can keep occurring repeatedly to producesmaller and denser objects until they are too dense tohave efficient radiative cooling. The mass at this “opacitylimit” is expected to be around m opacity ∼ M Jupiter . Thusthe opacity limit mass, rather than the initial core Jeansmass, would set the mass of the smallest fragments. Inthis scenario, it normally assumed that these small massesform first. These so-called embryos then accrete massand ultimately form both sub-stellar and stellar objectsin a scenario referred to as “competitive accretion”. Theaccretion of the surrounding material onto embryos isassumed to be rapid so that embryos remaining in denseparts of the molecular cloud would ultimately becomestellar-mass objects. Bate et al. (2003) explored this sce-nario with a simulation of cluster formation from a 50 M ⊙ molecular cloud and found embryos could remain inthe cloud undergoing rapid accretion or get ejected duringviolent interactions, resulting in a large number of low massobjects. More recent work indicates that a more realistictreatment of radiative energy losses can strongly suppressthe formation of sub-stellar objects in such simulations(Price & Bate 2009).The existence of extended disks around young stellarobjects is supported by both observations and theory (e.g.Yorke & Bodenheimer 1999). The disks are initially largeand flared at the class 0-I stage (extending to several thou-sands AU (see, e.g., Mundy et al. 2000) and are observableat mm or sub-mm wavelengths. They evolve into smaller,thin proto-planetary disks at the class III stage that aredetectable with infrared observations. Disk lifetimes are in-ferred to be about 10 Myr. Disks provide a natural environ-ment for forming sub-stellar objects because the Jeans massin disks is much smaller than in cores. It may be possible toform brown dwarf mass companions through the fragmen-tation of marginally stable ( Q ∼ Q > . M ⊙ and a half-mass radius, R . = 0 . × years (the approximate life-time for an early extended disk) with an impact param- eter at pericenter of 500 AU is about 30 percent. If sev-eral sub-stellar objects form during each encounter, as inWatkins et al. (1998a), a substantial numbers of BDs canbe produced. In this work, we focus on the outcomes of spe-cific encounters rather than estimates of their likelihood.The prior simulation work failed to address severalimportant issues. Firstly, the resolution used in the workcited above did not ensure that the Jeans mass was re-solved in the original disks. Thus it was not clear whetherthe objects might have condensed due to numerical effects(Bate & Burkert 1997; Truelove et al. 1997).Secondly, the range of circumstances in which gravi-tational instabilities can be triggered in a stable disk by anencounter were not explored. The Toomre criterion for insta-bility requires that that Q < Q crit for thin disks to be unsta-ble to axisymmetric perturbations, where Q is the Toomreparameter, locally defined according to (Toomre 1964), Q ≡ c s κπG Σ , (1)where c s is the sound speed, κ is the epicyclic frequency(close to the angular speed Ω for a Keplerian disk) and Σis the disk surface density. Analytical studies by Toomre(1964) and numerical studies of relatively thin proto-planetary disks (Boss 2001; Mayer et al. 2004; Pickett et al.2003; Gammie 2001) have shown that Q crit is in the range1-1.5 for various modes of fragmentation. If proto-stellardisks are stable according to this criterion then they tendto stay so indefinitely in the absence of mechanisms to re-distribute material. A non-negligible (vertical) thickness forthe disks can be stabilizing, reducing the critical Q valuefor axisymmetric perturbations from 1 down to 0 . − . c (cid:13) , 000–000 ncounter-induced brown dwarf formation cal simulations. This allowed the Jeans mass to be resolveduntil well after fragmentation initially took place and thevertical structure of the initial disks was also well resolved.The paper is organized as follows: in section 2 we describethe numerical method, the construction of the initial diskmodels and a study of resolution requirements. In section 3we present the encounter parameter space being explored,including Q values, relative orientations of disks and varia-tions in encounter velocity. Section 4 describes several typ-ical simulations and discusses the conditions for encounter-induced fragmentation. Section 5 presents analysis of theproperties of the resultant objects, including shapes, angu-lar momenta and orbital evolution. Section 7 summarizesthe results and includes discussion of the observational im-plications and possible directions for future work. Observing the exact mass, density structure and kinematicsof proto-stellar objects during the earliest stages is difficult(White et al. 2007). Using a hydrodynamical calculation ofproto-stellar disk formation, Yorke & Bodenheimer (1999)found that early disks contain a large fraction of the totalmass of the system and collapse relatively slowly. Based onthese results, we assumed that the proto-stellar disks werequasi-hydrostatic in order to construct a self-consistent diskdensity structure.We chose to simulate systems in the solar mass rangewith a single model for the disk-proto-star system where thedisk mass was about 0.6 M ⊙ , with a stellar mass of around0.5 M ⊙ . The disk extended to about 1000 AU. This aspect ofthe set-up is similar to the disks simulated by Watkins et al.(1998a,b). A power-law surface density profile Σ = Σ − p waschosen with p = 1 . max ∼ . g/cm .Within 100 AU the gas was initially smoothly truncated. Werelaxed our initial disks in isolation for about 10000 years,during which a small amount of material occupied the trun-cated region. Beyond 500 AU the surface density falls rapidlybut smoothly, and at 1000 AU Σ ∼ . g/cm . The entireprofile after relaxation is similar to the results from core col-lapse calculations by Yorke & Bodenheimer (1999). Most ofthe disk mass (0.5 M ⊙ ) is within 500 AU.It was assumed that the central proto-star was themajor heating source for the disk, and that the disk re-radiates like a black-body, giving a disk temperature pro-file of T ( r ) ∝ r − / (Mayer et al. 2004; Pickett et al. 2003;Chiang & Goldreich 1997) for the majority of the simula-tions. The profile smoothly transitions to an external, ambi-ent temperature of ∼ −
80 K. This profile is similar to theYorke & Bodenheimer (1999) results. With the radial sur-face density and temperature profiles specified, the verticalstructure was iterated to a self-consistent hydrostatic equi-librium taking into account both the gravity of the centralstar and the self-gravity of the disk. Accurate treatment ofthe gas self-gravity was important given the large disk massrelative to the star.
The choice of radial profile gave a nearly constant Q valuebetween 100-500 AU, with a minimum value ( Q min ) rangingfrom 1.6 to 2.1 depending on the case studied. The Q valuerises in the inner 100 AU and the disk is very stable there.Test simulations of single disk evolution with different val-ues of Q min were performed, and it was found that for ourdisk model values of Q min below Q crit = 0 . Q crit = 1 . Q min values used were at least twice as the Q crit for an isolated disk, i.e. Q min > . A fixed temperature profile T ( r ) ∝ r − / for most of thedisk was used in the majority of our simulations, whichis often referred to as a locally isothermal equation ofstate (EOS). Heating sources other than the star, such asshocks, were not taken into account. For the current study,it was found that the isothermal EOS was a good approx-imation for the initial disk and also acceptable even afterdisk was heated by tidal structures and shocks (outside theclumps) because the disks are extended and initially opti-cally thin. Given the temperature of the gas, using the sur-face density profile described above and adopting the opac-ities from D’Alessio et al. (2001), the product of the diskcooling time t cool and angular rotation rate Ω was calcu-lated and Ω t cool = 0 . − . t cool <
12 (Rice et al. 2005;Johnson & Gammie 2003) before the disk is dense enoughfor the isothermal assumption breaks down. However, afterthe disk fragments and the resultant objects became grav-itationally bound, the final densities of these objects weretypically orders of magnitude larger than the initial valuesand the objects were thus optically thick. Thus our locallyisothermal approximation was initially good but later be-came an upper bound on how efficient cooling should be.To further examine the importance of cooling in fragmenta-tion, a case was also investigated with explicit cooling thatbecame negligible in the optically thick regime. This secondassumption is a lower bound for the cooling in a real diskas it disallows even the slow leakage of radiation from densegas. The results of this case are described in section 6.
We used the TreeSPH code GASOLINE (Wadsley et al.2004) to model the disks during encounters. Around 200,000particles were used for each 0 . M ⊙ disk, with each particlerepresenting ∼ N neighbor (Bate & Burkert 1997)to prevent numerical fragmentation, where N neighbor is thenumber of neighbor particles in SPH (Smoothed Particle Hy-drodynamics) simulation. The number of particles ensured c (cid:13) , 000–000 S. Shen et al. that the Jeans mass was not-only well resolved in the initialcondition, but also well-resolved when density was enhancedduring an encounter (when the Jeans mass decreased), upto a density ρ ∼ − g cm − (about 10 times the ini-tial maximum value in the disk). Resolving the Jeans masswithin the fragments, however, was not attempted due tohigh computing costs ( N > particles usually required).Therefore, our choice of particle numbers ensures physicallycorrect modeling of the initial fragmentation in the proto-stellar disk, but cannot ensure that any secondary fragmen-tation of the clumps themselves are physical.The initial disk model also resolved the scale height atthe mid-plane very well at 100 AU (where the local smooth-ing length was about 1.9 AU and the scale height was about6 AU). Within 50 AU the scale height of the disk was lesswell resolved and the orbital time was shortest so that artifi-cial viscous evolution could have been a concern. The artifi-cial viscosity was moderated using the Balsara switch whichhelps to suppress unwanted shear effects. The measured vis-cous time scale in the code was about 10 to 10 years, muchlonger than the dynamical scale of the encounter, t ∼ years. Hence, the effect due to the viscous evolution was neg-ligible. Appendix A presents as exploration of the effects ofunresolved Jeans mass and unresolved disk scale height.The gravitational softening, ǫ , for gas particles was 0.2AU, one tenth of the smallest particle spacing in the initialcondition. Thus dense regions up to 1000 times the initialdensest part were well resolved. Since a larger ǫ usually sup-presses gravitation forces, the choice of gravitational soft-ening will tend to inhibit fragmentation once the particleseparations are comparable to or smaller than the soften-ing. In all tests, including those done for this work and inother work to date, there have been no indications that theGasoline code suffers from artificial fragmentation as a con-sequence of the gravitational softening being smaller thanthe local SPH smoothing length. The fixed softenings en-sure accurate gravitational dynamics at the cost of havingsomewhat conservative time steps in lower density regions.The central proto-star was represented by a single 0.5solar mass sink particle for which the softening was set to 1AU. The sink acted within a 1 AU radius around the starto consume gas falling into the innermost regions and thusavoid the small time steps required for gas orbiting veryclose to the star. However, no objects formed during thesimulation were replaced with sink particles, ensuring thatany disks surrounding the newly formed, sub-stellar objectwere modeled directly. In each pairwise disk encounter, the two disks were identicalin size, mass and Toomre parameter with an initial impactparameter of 1000 AU. Due to gravitational focusing duringrelatively slow encounters the minimum distance betweentwo stars could become substantially less than the initialimpact parameter. The initial disk separation was 6000 AU.This system was initially bound if the relative velocity ofthe encounter is less than 0.8 km/s and the final impactparameter is a function of the relative velocity. This studyexplored six dimensions in the overall parameter space: (i) The initial Toomre parameter profile characterized us-ing the minimum value, Q min ;(ii) The geometric configuration of the encounter, includ-ing the angles between the two disk planes and the anglesbetween the disk planes and the plane of the orbit of thetwo disks (a total of four independent parameters);(iii) The relative velocity of the encounter.Mass, disk size and impact parameter variations werenot included to limit the study to a manageable size. Smallerimpact parameters are possible but considerably less likely.These variations may be explored in future work.A list of all the simulations is shown in Table 1. Thefirst 3 cases (Q1, Q2, Q3) were carried out to study the effectof the Toomre parameter on encounter-triggered fragmenta-tion with different initial minimum Q values of 1.6, 1.8 and2.1, respectively. A higher Q profile means a more stableinitial condition, for which more effort is required to trig-ger fragmentation. We selected an encounter configurationwhere fragmentation was likely (coplanar, retrograde disks)to first explore the threshold in initial Toomre Q min , abovewhich dynamical interactions would become insufficient totrigger fragmentation. This study helped to constrain theregion of interest within the parameter space.Previous work by Lin et al. (1998) indicated that frag-mentation can take place in tidal structures forming duringencounters while Watkins et al. (1998a,b) suggested diskswould fragment. The degree of tidal distortion and angu-lar momentum transport has been shown to depend on theinitial disk spin orientations (with respect to the orbital an-gular momentum) (Toomre & Toomre 1972). Various diskspin orientations were used in our simulations including pro-grade cases (where the angle between the angular momen-tum vector of the disks and that of the orbit was less than90 degrees), retrograde cases (where the spin-orbit angle wasgreater than 90 degrees) and combinations of the two.In most simulations the disks were not close to the sameplane, as would be expected for randomly oriented disks.However, coplanar encounters (where the spin-orbit anglesare zero) were also simulated as special cases to investigatethe shocks generated in such encounters and the possibilityof shock fragmentation. In every case the orbital plane wasthe x-y plane so that the orbital angular momentum was inthe direction of the z-axis.In a completely generic encounter, as shown in Figure 1,the disks can be oriented randomly. The configuration ofeach disk in the encounter can then be characterized by twoangles: the angle between the spin angular momentum ofeach disk and the positive z direction, θ , and the angle be-tween the positive x direction and the projection of diskangular momentum onto the x-y plane, φ , which also deter-mines where that disk intersects the x-y plane. The initialdisk velocities are always along the x-axis and the initialimpact parameter is due to an offset in the y-direction.Cases Conf1-48 explored the ( θ , φ , θ , φ ) parameterspace. For a single disk, the angle parameters were selectedsuch that the θ space (from 0 to 180 degrees) was dividedevenly into three ranges and the φ space (from 0 to 360degrees) into four for a total of 12 cases. Naively there wouldhave been 144 possible combinations of disks but due to thesymmetry of the encounters and the equivalence when thetwo disks were exchanged in space, there were 42 unique c (cid:13) , 000–000 ncounter-induced brown dwarf formation xz φ y +v/6000 AU1000 AU −v/ φ θ θ Figure 1.
A generic encounter configuration between two proto-stellar disks. The centres of mass orbit in the x-y plane with an initialvelocity solely in x-direction. The angular momentum of an individual disk has a vector direction that differs from the orbital angularmomentum (the z-axis) by an angle θ . The disk angular momentum vector projected onto the x-y plane makes an angle φ with the x-axis. collisions to simulate. An additional six cases were simulatedwhere at least one disk was exactly in the orbital plane, toinvestigate the effects of coplanar encounters.The cases with labels starting with “Vel” refer to twoseries of runs used to examine the constraints due to the ini-tial relative velocity on disk fragmentation for prograde andretrograde encounters. Relative encounter velocities in therange of 0.8 km s − to 6.0 km s − were used in the simula-tions to investigate the effects of relative velocity on diskfragmentation. The inferred disk velocities are consistentwith the typical velocity dispersion in stellar clusters (a few km s − , See e.g., Bate et al. 2003, and references therein) To explore the impact of the initial minimum Toomre pa-rameter we investigated 3 cases, all of which have ( θ, φ ) =(180 ,
0) and relative v = 0 . θ and φ .Case “Q1” (Q min =1.6) is depicted in Figure 2. Whenthe physical impact began, the gas was shock compressedand quickly formed a “shock layer” between the two diskswhere the density was substantially enhanced. At 24,000years after the encounter, the shock fragmented to form a7 M J clump at about 223 AU away from the nearest star(label “a” in Panel 5, Panel 7). Also, due to its retrogradeconfiguration, the angular momentum of the disk spin andthe orbit are opposed and thus the rotation speed of thegas decreases. As a result, a large amount of gas spirals intothe inner disk (within ∼
100 AU in this case), while gas in the outer disk dissipates. The azimuthally averaged surfacedensity profile of the lower-left disk in Panel 4 of Figure 2 isshown at different times in Figure 3 and the correspondingevolution of azimuthally averaged Q profile is shown in Fig-ure 4. At t = 23,000 years, the surface density had increasedsignificantly, and at r ∼ − AU , it had almost doubledits value from earlier times (t = 6,000 years, well before thephysical impact of gas materials). However, Q was still veryhigh in the inner 100 AU because of the fixed high temper-ature in that region (Figure 4). Thus it was the presenceof the shock layer at around 200 AU that locally increasedthe surface density and lowered the Toomre Q down to unity( Q min = 1, Figure 4) allowing a strong spiral feature to formwhich is gravitationally unstable. At t = 24,000 years, a pairof close clumps formed at r ∼
130 AU with masses of about12 and 14 M J (labeled “b” in Panel 5, Panel 8), and atabout 150 AU another gas blob started condensing (a mod-erate extremum near the 130 AU clump in the plot of the Qprofiles, labeled “c” in Panel 5), which eventually became abound object (labeled “d” in Panel 6). Gas inflow increasedafter formation of these objects. At t = 24,400 years, frag-mentation within 20 AU was finally triggered in one of thedisks and 4 brown-dwarf mass clumps were formed (boxedregion in Panel 6, Panel 9). However, these objects were notincluded in our statistics (Section 5) because the densitywithin 20 AU at this time step exceeded 10 − g cm − andbecame optically thick (cf. Section 6), so that an isothermalcooling assumption was no longer appropriate. The simula-tion was stopped at t = 24,400 years because the time stepreduced dramatically after the first condensation formed.Analysis based on the energy and momentum at the end ofsimulation indicates that amongst the 4 objects formed out-side 100 AU , one is unbound and likely to leave the system,one has a large semi-major axis (a = 879.3 AU) and theother two have orbits with semi-major axes less than 100AU. The encounter “Q2” has same parameters as “Q1” ex-cept the initial Q min = 1 .
8. The disk underwent similar c (cid:13) , 000–000 S. Shen et al.
Table 1.
Simulation ParametersCase Toomre Q min ( θ , φ ) (deg) ( θ , φ ) (deg) v (km/s) N N Q1 (Vel retro1) 1.6 (180 ,
0) (180 ,
0) 0.8 3 (1) 0 (4)Q2 1.8 (180 ,
0) (180 ,
0) 0.8 0 0 (1)Q3 2.1 (180 ,
0) (180 ,
0) 0.8 0 0Conf1 1.6 (180 ,
0) (180 ,
0) 2.0 0 (1) 0 (1)Conf2 1.6 (0 ,
0) (0 ,
0) 2.0 3 4Conf3 (Vel pro2) 1.6 (0 ,
0) (45 ,
0) 2.0 1 2Conf4 1.6 (180 ,
0) (135 ,
0) 2.0 1 0 (1)Conf5 1.6 (180 ,
0) (45 ,
0) 2.0 2 (1) 0Conf6 1.6 (150 , , , , , , ,
70) (150 ,
70) 2.0 2 3Conf10 1.6 (150 ,
70) (150 , ,
70) (150 , ,
70) (150 , , , , , , , ,
70) (150 ,
70) 2.0 0 3 (2)Conf17 1.6 (30 ,
70) (150 , ,
70) (150 , ,
70) (150 , , , , , , , ,
70) (30 ,
70) 2.0 1 (1) 1 (1)Conf24 1.6 (30 ,
70) (30 , ,
70) (30 , ,
70) (30 , , , , , , , ,
70) (90 ,
45) 2.0 0 3 (1)Conf31 1.6 (30 ,
70) (90 , ,
70) (90 , ,
70) (90 , , , , , , , ,
45) (150 ,
70) 2.0 0 2 (1)Conf38 1.6 (90 ,
45) (150 , ,
45) (150 , ,
45) (150 , , , , , , , ,
45) (90 ,
45) 2.0 0 0Conf45 1.6 (90 ,
45) (90 , ,
45) (90 , ,
45) (90 , ,
0) (90 ,
0) 2.0 0 3 (1)Vel pro1 1.6 (0 ,
0) (45 ,
0) 0.8 0 (1) 0Vel pro2(Conf3) 1.6 (0 ,
0) (45 ,
0) 2.0 1 2Vel pro3 1.6 (0 ,
0) (45 ,
0) 4.0 0 0Vel pro4 1.6 (0 ,
0) (45 ,
0) 6.0 0 0Vel retro1(Q1) 1.6 (180 ,
0) (180 ,
0) 0.8 3 (1) 0 (4)Vel retro2 1.6 (180 ,
0) (180 ,
0) 2.0 0 (1) 0 (1)Vel retro3 1.6 (180 ,
0) (180 ,
0) 4.0 1 (1) 0Vel retro4 1.6 (180 ,
0) (180 ,
0) 6.0 0 0The pairs ( θ , φ ) and ( θ , φ ) denote the angles between the disk angular momentum vectorand the angular momentum vector of the orbit and the angular momentum vector and thevector direction of the initial velocity difference in the orbital plane for disk 1 and disk 2respectively. v denotes the relative encounter velocity. The number of objects formed foreach disk are shown in the last two columns. If objects were formed within 50 AU then thenumbers of these are shown in parentheses. c (cid:13) , 000–000 ncounter-induced brown dwarf formation −13 −6 −13 −6 −13 −6−10 −4 −10 −4−13 −6−10 −4 −10 −4 −9 −3t = 6,000 yrs t = 12,000 yrs t = 18,000 yrst = 24,400 yrst = 24,000 yrst = 23,000 yrs 1 2 34 5 6987 bac d Figure 2.
A coplanar encounter between two retrograde disks with Q min = 1.6 . Panels 1 to 6 are snapshots of the simulation, withthe time indicated at the left top of each panel. The disk were rotating clockwise in all of the panels. The color bar under each panelgives the logarithm (base 10) of the density in units of M ⊙ /AU . In panels 4 and 5 only the inner dense regions were included and adifferent density range was used as indicated. Panels 7 and 8 are zoomed-in pictures of objects formed from a shock layer and due to diskfragmentation (a and b in panel 5), respectively. Panel 9 is the zoomed-in picture of the boxed region in panel 6, showing fragmentationof the inner disk.c (cid:13) , 000–000 S. Shen et al.
Figure 3.
Evolution of the surface density profile for the case inFigure 2.
Solid line : Surface density profile of one of the encoun-tering disks at 6,000 years.
Dashed line : Surface density profilefor the same disk at 23,000 years.
Dotted line : Surface densityprofile at 24,000 years. The sharp peak at 130 AU is produced bytwo close objects forming due to disk fragmentation, the one at225 AU originates from shock layer fragmentation. The moderatepeak in between indicates ongoing gas condensation.
Figure 4.
The Toomre Q profile for the same disk as in Figure3 at 6,000, 23,000 and 24,000 years (solid, dashed, dotted line,respectively). The dot-dashed line indicated Q = 1.0 physical processes such as formation of a layer of shock-compressed gas, inflow of gas and disk truncation after theencounter. Nevertheless, the disks were more stable to theperturbation and only one of the disks developed an m = 2mode gravitational instability at a relatively late time step(t = 25,300 years), which then fragmented and formed onlyone bound object at radius of 30 AU. The other disk wastruncated but remained stable. The gas in the shock layerdid not condense to form an object. Instead, it dispersed af-ter the disks moved apart. The lowest Q value at t = 24600years (the time step just before fragmentation) for the un-stable disk was 1.1 and for the stable disk was 1.7, close tothe original minimum.The last case in this series “Q3” used initial disks with Q min = 2.1. The encounter disrupted both disks withoutforming any objects. The two stars happened to capture each other and formed a close binary system with a separation ofabout 20 AU. Most of gas was dispersed, with the bound gasat the end of simulation forming a circumbinary disk witha radius of 250 AU at 31000 years. The circumbinary diskhad a much lower surface density compared to the initialdisks. Assuming the gas was still in Keplerian rotation, theminimum Toomre Q was very high ( Q min = 6 . ∼ K ). Sogas inflow can be the key mechanism in encounter-induceddisk fragmentation.Note that the local enhancement of density in theshocked gas can lower the Q parameter significantly andinduce disk instability before the inflowing gas reduces Qin the inner, hotter region, making the fragments form rel-atively far from the stars. This is especially clear for copla-nar encounters discussed in this section. In more likely, non-coplanar cases (cf. Section 4.2) the shocked gas usually doesnot form a well-defined layer between the disks, but it re-mains an important factor driving disk instability. In thecase “Q1” the shock layer itself also fragmented directly.The object formed this way was not associated with spiralstructure in the unstable disk (panel 5 in Figure 2). Frag-mentation of shock layers rarely occurs and was restricted tocoplanar, low-velocity, low-Q cases. With higher disk tem-perature (as in “Q2”) or higher velocity the shock layer usu-ally dissipates when the two stars pass periastron withoutforming objects. The parameter space of encounter configurations was ex-plored with 48 encounters covering the ( θ, φ ) parameterspace. The effect of the relation between disks’ spin and or-bital angular momentum in z direction can be see by varyingthe disks’ θ parameters. This variation qualitatively changesthe interactions, including the mechanisms of fragmentation.There is a large stochastic element to the evolution ofthe disks and the formation of clumps. In this sense, the pre-cise number of objects formed in a single simulation shouldnot be examined too closely. For example, a small perturba-tion in the initial particle placements that is within the noisein the initial surface density profile is sufficient to change theoutcome in terms of the precise number and placements ofclumps. For example, several cases are symmetric from thepoint of view of the two disks but have non-symmetric out-comes. Symmetry in the configurations is revealed by themappings x, y, z → − x, − y, z which implies θ = θ and c (cid:13) , 000–000 ncounter-induced brown dwarf formation φ = φ (such as for the coplanar cases and “Conf6”, forexample) and x, y, z → − x, − y, z which gives θ = θ and φ = φ + 180 (e.g. “Conf11”). These initial states werecreated using absolute offsets ( x → x + ∆ x ) of two iden-tical disks that have small non-axisymmetric perturbationsdue to the glass initial particle state. Thus the small pertur-bations in the particle configurations are not symmetric inthese cases. These difference magnify during the course ofthe encounter and ultimately result in different numbers ofobjects in each disk as can be seen in the object numbers foreach disk in symmetric cases in table 1. For example disk 1in “Conf6” made 2 objects while Disk 2 made 3. The vari-ations are consistent with a Poisson-type random process.In most of our analyses we look at outcomes averaged overseveral encounters. Encounters between retrograde disks tended to reduce thedisk spins and cause gas inflows, increasing the surface den-sity of the inner disks and reducing the Toomre Q parameter(Figure 2). Prograde encounters are more complex as seenin Figure 5 showing snapshots of encounter “Conf3”. Thetwo disks have ( θ, φ ) = (0 ,
0) and (45 , ∼ M J .A tidal structure was also present in the other disk, wherethere was a 45 degree angle between the spin and orbitalangular momentum vectors. At 15,000 years, it condensedinto two objects with masses of 9 M J and 70 M J at 60 and110 AU from the star, respectively (labeled “b” in Panel 5Figure 5). The simulation was terminated at 16,000 years.Analysis based on the energy and momentum of the objectsfound that all the 3 objects remained bound to the systembut their orbits were large: 532 AU for the one that con-densed first, and 90 and 307 AU for the two condensed inthe later time step.Fragmentation of tidal structures was seen in severalprograde disk encounters, especially when the disk wasclosely aligned with the orbital plane. This result is consis-tent with Lin et al. (1998), in which a disk-disk encountertriggered formation of a large tidal tail and an object con-densed from it. Objects formed this way tend to have largeorbital radii. For example, in the case“Conf2”, tidal struc-tures from both disks condensed and formed 7 objects intotal by the end of simulation, all of which were unbound toboth stars. Figure 6.
The number of objects formed in a disk per encounteras a function of θ . Each box corresponds to a single simulated diskthat produced the number of objects indicated on the y-axis. Theboxes have been spaced apart in the x-direction for clarity. Theaverage and variance are indicated with solid and dashed linesrespectively. Tidal structures due to spin-orbit resonances are not as ef-ficient at producing objects as the inner disk fragmentationseen in retrograde disks. On the other hand, Watkins et al.(1998a) has argued that star formation triggered by col-lisions of molecular clouds tends to result in prograde-prograde coplanar encounters so they may be more common.In Figure 6 we plot the number of object formed in adisk per encounter as a function of θ , the angle between theindividual disk spin vector and the positive z direction. Wechose the homogeneous sample from “Conf6” to “Conf48” inTable 1, excluding the unlikely coplanar encounters, φ = 0encounters and encounters with disks residing in the orbitplane. The cases chosen all have same initial Toomre param-eter and relative encounter speed.The average number of bound and unbound objectsformed within that disk increases with increasing θ for thatdisk. This is the case regardless of the configuration of theother disk in the interaction. Gas inflow efficiently adds massto inner disks which fragment whereas in prograde disks, themass transfer usually send substantial amounts of gas out-ward rather than inward. The condensation of tidal struc-tures did not occur frequently, because in most of the casesthey usually dissipated before anything could form. No clear correlation was found between the productivity ofa disk and the disk’s own φ parameter. However, no frag-mentation would happen if the two disk planes are parallel(or nearly parallel) initially but not coplanar. 60% of theencounters in the “Conf” series which do not produce anyclumps fall into this category. The disks were simply trun-cated in these cases. The final size of the disk is about 300to 500 AU depending on the initial configuration. c (cid:13) , 000–000 S. Shen et al. t = 12,000 yearst = 10,000 years t = 13,000 yearst = 14,000 years t = 15,000 years t = 16,000 years−12 −5 −12 −5 −11 −4−6−6 −13−6−13−13 a b
Figure 5.
Snapshots of the prograde-prograde disk encounter “Conf3”. Time is indicated on the left top of each panel. The disks arecounter-clockwise rotating in this view. The color bar under each panel gives the logarithm (based 10) of the density in units of M ⊙ /AU .Different density ranges are used to clearly show the clumps. We investigated the effect of varying the relative velocitywith two series of simulations, one containing non-coplanarencounters with both disks prograde, and one with coplanarencounters with two retrograde disks. In each series, the rel-ative velocity of the encounters is repeated along with thenumber of object formed in Table 2.With the exception of a special single case (“Vel pro1”,see below), the number of clumps formed in an encountergenerally decreases with increasing velocity. No clumpsformed in prograde disks when the relative encounter ve-locity was above 4.0 km/s. This result can be qualitativelyunderstood by comparing the interaction timescale with thedynamical timescale of the disks. Though gravitational in-teractions occur regardless whether or not there is a physicalimpact, these simulations indicate that significant changesin the density profile occur only when the disks are closeenough to physically impact. The interaction timescale is
Table 2.
Number of objects formed in the simulations with dif-ferent encountering velocitiesCase Velocity Number Fragmentation(km/s) of Clumps MechanismVel pro1 0.8 1 Circumbinary diskVel pro2 2.0 3 Tidal structure and diskVel pro3 4.0 0Vel pro4 6.0 0Vel retro1 0.8 8 Shock layer and diskVel retro2 2.0 2 Disk onlyVel retro3 4.0 3 Disk onlyVel retro4 6.0 1 approximately, t int = d eff v (2)where d eff is the twice the disk diameter (4000 AU in this c (cid:13) , 000–000 ncounter-induced brown dwarf formation case), which gives the rough separation for physical impact. v is the relative velocity of the encounter. The dynamicaltimescale within each disk is a function of radius r , givenby, t dyn ( r ) = 2 πr / √ GM ∗ (3)where G is the gravitational constant and M ∗ is the massof the star. The encounter timescales t int were 2 . × ,9 . × , 4 . × and 3 . × years for relative speeds 0.8,2.0, 4.0 and 6.0 km/s, respectively. The dynamical timescale t dyn was about 1 . × years at 400 AU and 3 . × yearat 200 AU. The implication is that there was insufficienttime during fast encounters, where t int < t dyn to apply thetorque that transfers angular momentum to change the diskprofiles and causes fragmentation. This is more restrictivefor prograde disks where tidal structures form further outat distances similar to 400 AU. In the case “Vel pro1”, star capture occurred and the twoencountering stars formed a binary system. Both disks wereprograde and large tidal tail structures formed when thedisks passed periastron. The tidal structure dissipated overtime and carried away the excess angular momentum fromthe system, enhancing the formation of a close binary. Theseparation of the two components was about 90 AU, signif-icantly smaller than the periastron of the trajectory (about300 AU taking into account the effect of gravitational focus-ing). The circumbinary disk was about 500 AU in radius.At a later time it also fragmented to form a 20 Jupiter massobject.Another case where stellar capture occurred is the case“Q3” (as described in Section 4.1). Since the starting disk Qparameter is 2.1, the circumbinary disk became hotter andless massive compared to the standard cases. No fragmentsformed at the end of the simulation. The initial disks areretrograde and the encounter did not induce tidal structures.The two cases above are the only ones where stellar cap-ture occurred. Note that both of them have rather conser-vative choices of encounter velocities (0.8 km/s), with whichthe system are marginally bound at the initial state.
Most of the encounters produced clumps with the excep-tion of those with high relative velocities or very stable ini-tial disks (e.g. Q min > Figure 7.
The mass distribution of 146 objects formed outside50 AU in the 48 “Conf” simulations. The bin size is 5 M J . the initial state, the disks were well resolved outside 50 AU.However, as the simulations progressed, additional matteraccumulated in the inner disks improving the resolution el-ements per scale height. On the other hand, objects with50 AU are likely to experience further evolution such as ac-cretion and migration into the star. In addition, work byRafikov (2007) has suggested that disks this close to a starare optically thick and unlikely to cool fast enough to frag-ment. For these reasons, we have excluded objects within50 AU from the sample discussed below. However, the ob-jects formed within 50 AU are not statistically different andproperties such as the relative mass distribution would notbe significantly altered if these objects were included. A to-tal of 191 objects formed in our simulations (excluding thehigh Q cases) with 144 outside 50 AU. The object masses range from 0.9 Jupiter masses ( M J ) to127 M J , extending into the stellar regime. The mass distri-bution is plotted in Figure 7. It was found that the numberof objects generally decreases with increasing mass, con-sistent with the observed sub-stellar initial mass function(IMF). The simulated population does not include objectsmuch below 10 M J ). In this regime there are 1000 or fewerparticles per object. The numerical gravitational softening(0.2 AU) could also inhibit the formation of smaller objects.Thus these simulations are unable to probe whether there isa physical mass cut-off in this range.Observationally, the sub-stellar IMF in open clusterscan be well described by a power-law dN/dM ∝ M − α .For young open clusters such as Pleiades, λ Orionis and σ Orionis, the exponent α is close to 0.6 at low masseswith an uncertainty of order 0.1 (Moraux et al. 2003;Barrado y Navascu´es et al. 2004; Caballero et al. 2007). Anuniversal IMF was proposed by (Kroupa 2001) in which α = 0 . ± . . M ⊙ < M ∗ < . M ⊙ . Inthis work, we fit the mass distribution of the sub-stellar re-gion ( M < . M ⊙ ) and the best fit gives α = 0 . χ uncertainty 0.15 (Figure 8), consistent with the observedvalues. Above the stellar boundary the number of objectsdeclines more rapidly with mass and a break is apparent c (cid:13) , 000–000 S. Shen et al.
Figure 8.
The logarithmic mass distribution of the objects. Theerror bars were calculated assuming Poisson errors ∝ √ N bin . Stars : sub-stellar objects (
M < = 80 M J ); Diamonds : low mass stars(
M > M J ); Dashed line : linear fit of the logarithmic distribu-tion of sub-stellar objects dN/dM ∝ M − α , with α = 0 . near 0.08 M ⊙ (or about 0.1 of the central mass) which isalso consistent with the trend reported by (Kroupa 2001).There are only 9 objects in total that are beyond 0.08 M ⊙ and thus the statistics are poor. More simulations, includ-ing those with varying central star masses, are required toexplore the resulting mass function in the low mass, stellarregime. Most of resultant clumps in our simulation are not sphericalbut have highly flattened, disky shapes (height to radius ra-tio ∼ to 10 g cm /s.Assuming the total angular momentum is conserved duringthe evolution and adopting a typical radius of a brown dwarfto be 7 . × cm (Burrows et al. 1998), the rotation speedat the edge of the brown dwarf would be 10 times the breakup speed. Thus the proto-brown dwarfs disks must shed ro-tation in order to collapse down to a proto-BD. This result isconsistent with recent mid-IR, sub-millimeter, and millime-ter observations of the spectral energy distributions (SEDs)of young brown dwarfs, which indicate that most of themare surrounded by circumstellar disks (e.g., Muench et al.2001). A T Tauri accretion phase has also been indicated insome young BD disks by the detection of broad, asymmet-rical Hα lines (Jayawardhana et al. 2003). The IR excess inSEDs was found to decrease as the age of BDs increases,which suggests that the BD disks evolve from flared to flatgeometry similar to proto-stellar disks (Mohanty 2005). Theclumps in our simulations are flat partly due to the assumedisothermal EOS. As explored in section 6, these dense con-densations are probably optically thick and would thus behotter and more puffed up in nature.The simulations were terminated a few thousand years −7 −2−7 −2 Figure 9.
A resultant clump with typical disky shape. Both pan-els are 14 AU on a side. The clump has disk radius about 4 AU.
Left : The object projected in the x-y plane;
Right : The objectprojected in x-z plane. The gray scale indicates the density inlogarithmic space in code units M ⊙ /AU − . after first object formed due to computer time limitations.In some objects, like the one shown in Figure 9, by the timethe simulation was halted, the gas with lower specific an-gular momentum had started to condense into the centerand a density gradient had also been established. Becausethe gravitational softening (0.2 AU) is comparable to thesize of most clumps, our ability to resolve the structure andfollow the evolution of the proto-BD disks is limited. How-ever, it can be expected in later evolution a central protoBD will form and the surrounding disk will accrete onto iton a viscous timescale.In the majority of cases where clumps formed fromfragmentation of disks, the spin directions of the proto-BDdisks are similar to their parent proto-stellar disks. In thecase where the shock layer fragments, however, the resultantproto-BD disk may spin in a direction different from bothdisks. For the proto-BD disks that are still orbiting the star,further interaction with their parent disk will influence thefinal spin. Typically where fragmentation occurred more than one ob-ject was produced. About 15 percent of the objects are inBD-BD binary systems. This fraction is consistent with theobserved binary fraction (around 10 ±
10 % at the range0.01 to 0.1 M ⊙ ) (cf. Figure 3 Hubber & Whitworth 2005).In most cases, the binaries formed from secondary fragmen-tation of very large disky clumps. Further numerical stud-ies with higher resolution are necessary to investigate thisquestion in detail because, although the Jeans mass waswell resolved before fragmentation and most of the clumpsare resolved by at least 1000 particles, the Jeans massesof the clump themselves decrease so quickly with increas-ing density that they are only marginally resolved in oursimulations. Without resolving the Jeans mass, it has beensuggested that fragmentation could be artificially enhanced(Bate & Burkert 1997). However, given the very high spinrate of the clumps, secondary fragmentation is likely to oc-cur, according to Matsumoto & Hanawa (2003).Most of the BD binary systems were still bound to thecentral stars when the simulations were terminated. For ex- c (cid:13) , 000–000 ncounter-induced brown dwarf formation ample, in case “Conf5”, the binary system orbits the starat ∼
140 AU, with the two components (mass 19 M J and13 M J ) orbiting each other in a period of around 30 yearsand with separation about 3 AU. Each component has ahighly flattened shape. In future evolution, it is expectedthat some of these binaries will merge to the star, leave thesystem, or be ejected in close interactions with other objects.The surviving systems are likely to remain hierarchical star-BD multiples. Thus the collision mechanism may explainthe origin of observed star-BD multiple systems such as GL569, where a recently confirmed BD triple (GL 569B) or-bits the primary, an M2.5V star (GL 569 A) (Gorlova et al.2003; Simon et al. 2006). The system is about 100-125 Myrsold so no BD disks were observed. However, in younger BDmultiples accretion signals have been detected (Kraus et al.2006), indicating accretion disks surrounding components ofthe multiple. Objects in our simulations typically formed at least 50 AUaway from the closest star. Approximately 60 percent of theobjects formed were not bound to the system. The percent-age is higher for coplanar encounters, in which the shocklayer fragmentation can take place relatively far out in thegravitational potential of the stars. The unbound objectsmay leave the whole system to become free floating BDs(or BD binaries) or low mass stars. For the objects that arebound to the star, the orbital properties vary from case tocase. About 5 percent of the orbits have very large semi-major axes ( >
500 AU) but are still bound to the system.Since the outer proto-stellar disk was dispersing in most en-counters, sub-stellar objects in large orbits were unlikely toaccrete additional gas to become low mass stars. Hence theseobjects may remain as wide-separation BD companions,which have been detected in observations (Gizis et al. 2001).Where more than one object orbits the star on nearby orbitsdynamical interactions can take place resulting in the pref-erential ejection of lower mass objects (Reipurth & Clarke2001). The 11 objects that were ejected to significant dis-tances during the simulations were lower in mass. However,the majority of the objects were unbound and would ulti-mately become free-floating.For the clumps that stay in the dense region of theproto-stellar disks (within 100 AU), the future growth andevolution depend on the rate that the gas accretes onto theobjects. Under the assumption of efficient gas accretion, as inthe competitive accretion scenario (Reipurth & Clarke 2001;Bate et al. 2003), these objects may eventually become low-mass stars over timescales much longer than those simulatedhere.
The simulated disks produced dense regions as they evolvedand reached the point where they were no longer uniformlyoptically thin. The disk material itself tended to remain op-tically thin ( τ <
1) and was thus capable of efficient cooling but the fragments within the disk reached high densitieswhere the radiative cooling times became very long. Thuson the timescale of these simulations, the interiors of thefragments were effectively adiabatic.In earlier work on proto-planetary disk fragmentation(Mayer et al. 2004), the approximation of converting the en-tire simulation to adiabatic gas was employed to test thebehaviour in this regime. This corresponds to the extremewhere cooling is completely absent. A less extreme test withstrong suppressed cooling is to smoothly interpolate fromefficient black-body cooling when the gas is optically thinto no cooling in the optically thick regime. Optically thinblack-body radiative losses are given by, L = 4 κσT M (4)where L is the luminosity, κ is the opacity, T is the tem-perature and M is the total mass. The opacity was takento be 1 cm g − at T = 60 K with a power-law variationwith temperature as T . . This is consistent with the dustopacity values calculated by d’Alessio et al., as tabulated inMejia (2004), over the range of 10-300 K.We then took a late stage of a representative run“Conf30” at 15000 years, after fragmentation had occurred,and measured the optical depth to infinity in 12 directions(corresponding the faces of a dodecahedron) for each parti-cle. This allowed us to estimate the averaged effective opticaldepth, τ eff = − log e R Ω e − τ (Ω) d Ω / π . We found that the ef-fective optical depth correlates fairly well with gas density inthe regime, ρ = 10 − to 10 − g cm − , where the materialtransitions from optically thin to ( τ eff ∼ .
1) to opticallythick ( τ eff ∼
10) as τ eff = 10 . . ρ . The rms deviationof the fit is 0 .
16 dec in τ . Assuming isotropic emission, wecan then estimate the fraction of the emitted light (givenby equation 4) that escapes the disk entirely. We then as-sume that the rest of the emission is absorbed on-the-spot,an approximation commonly used for optically thick lines inmodels of emission nebulae, so that the net radiative loss, L p , from one particle is L p = 4 κσT p m p e − τ eff ( ρ p ) , where m p , T p and ρ p are the particle mass, temperature and density,respectively. In the limit of infinitesimal particle masses, thisexpression is proportional to the contribution to the surfaceflux from a mass of gas, dm = ρ dl dA = dτ /κ dA . In the caseof an infinite plane-parallel slab, when an analytic expres-sion for τ eff is used and the flux contribution is integratedfrom the surface to high optical depth the total surface fluxasymptotically approaches the value of σT p .The assumption for the main set of simulations was thatthe disk temperatures were set by heating from the centralstar offset by radiative losses so that the net temperatureswere a function of the distance to the nearest star. To includethe effect of heating by the star, the governing equation ofthe particle energy per unit mass, E p , is, d E p dt = 4 κσe − τ eff ( ρ p ) ( T ( r star ) − T p ) − Pρ ∇ . v | p + Λ shock ,p , (5)where T ( r star ) is the temperature determined from an equi-librium between stellar heating and radiative losses at theradius r star from the nearest star, − Pρ ∇ . v is the compres-sive heating term and Λ shock is the shock heating term.Thus the model now includes compressive heating and shockheating. Outside the intense collision environment, wherethese process are unimportant in relation to stellar heat- c (cid:13) , 000–000 S. Shen et al. ing, the disk temperatures relax to the temperature asa function of radius assumed in the preceding sections.Overall, this similar to the approximation employed byStamatellos & Whitworth (2009), where local cooling wasbased on the amount of radiation escaping from regions withsimilar physical conditions in spherical collapse models. Inthe approximation used here the radiative loss rates werecalibrated from similar disk cases.Our approximation retains thermal energy indefinitelyin the denser gas and is thus conservative with regard toensuring that the fragmentation is not a result of favourablecooling assumptions. Detailed radiative transfer calculationsshould give results between this approximation and theisothermal EOS. It has been argued by Rafikov (2007) thatthe optically thick parts of proto-stellar disks are likely tobe convective, in which case our neglect of radiation as aheat transport mechanism may not have a large impact onthe disk structure. However, the cooling of the optically thinsurface layers are approximately correct in this model. Fullradiative transfer has been implemented within the Gaso-line code (Rogers, Shen & Wadsley, in preparation ) and willbe applied to proto-stellar disk collisions as part of futurework.
We re-ran the simulation we used to calibrate the effectiveoptical depth-density relationship using the new cooling ap-proximation. Figure 10 compares the disks at 14,000 yearsin the locally isothermal run (Panels that are labeled with“ISO 1” and “ISO 2”) with the ones in the new run (Pan-els labeled with “THICK 1” and “THICK 2”). Using theisothermal equation of state, one of the disks fragmentedand formed ten bound objects, including two binary sys-tems. However, when using the optical thick approximation,the cooling is suppressed and the disk does not fragment im-mediately. Instead, it is the tidal structure induced by thegravitational resonance that condensed and formed a clumpat about 410 AU from the star. At later time steps, the gasin the tidal tail dissipated and the objects orbit the star atabout 350 AU. The clump had a mass of 0.075 M ⊙ , closeto the hydrogen burning limit.Compared to the objects that formed in the correspond-ing isothermal run, the clump had a more spherical shape.Nevertheless, the angular momentum of the object was stillvery high, about 7 . × g cm /s. By the end of the simu-lation we did not see flattening as the clump could not coolwith this approximation. If no angular momentum was lostthe final object would spin at 30 times the break up speed.Hence we expect the clump to collapse to a disk similar tothe ones formed in the isothermal run and that jets andfragmentation could occur in nature as before.Thus with heavy cooling suppression, beyond that ex-pected in nature, we retain our main result, i.e. there wasstill gas condensation and the formation of clumps due toencounters. The number and the shape of the clumps atthe end of simulation was strongly affected. Exploring theeffect of differing thermal regulation on encounter-inducedfragmentation will be a target for future work. −11 −5 −11 −5−10 −4 −10 −4−10 −4 −10 −4ISO_1ISO_2 THICK_1THICK_2THICK_clumpISO_clump Figure 10.
A Q = 1.6 disk encounter with isothermal EOS vs.our optical thick approximation.
Left : The run with isothermalEOS (top and middle panels) and the typical shape of a clump(bottom panel);
Right : The simulation that used optical thickapproximation. The grey scale under each panel gives the densityranges in units of M ⊙ /AU . In this work, the possibility of encounter-induced disk frag-mentation as a mechanism for forming sub-stellar objectswas investigated. The significant increase in resolution overprevious work and the use of a realistic initial disk modelensures that the initial fragmentation to form clumps wasphysically modeled within the constraints of our thermalregulation assumptions. These simulations have shown thatencounters can trigger the formation of brown dwarfs as wellas low mass stars.We found that strong instabilities can be triggered instable disks with initial Q min < .
1, by lowering the local Qvalue with large gas inflows during an encounter. A relativelysteep surface density profile (Σ ∝ r − / ) and a relative highambient temperature (40 K) were used in our simulations. Ifdisks have shallower density profiles, as was found in someobservations of early disks (Andre et al. 2000), and hence c (cid:13) , 000–000 ncounter-induced brown dwarf formation have large portions of their mass in the cooler, outer region,the rate of fragmentation could be even higher.By studying collisions with random initial orientations,it was also found that retrograde encounters (disks with ret-rograde spins relative to the orbital angular momentum)more effectively transport angular momentum, triggers gasinflows and ultimately produce more clumps on average thanthe prograde ones (Figure 6). These results differs from thosereported for simulations of star-disk encounters which indi-cated that prograde encounters were more efficient in trig-gering disk instability (e.g. Boffin et al. 1998; Lodato et al.2007). The difference may relate to having an extended per-turber.Dense shock layers may occur during coplanar encoun-ters and result in objects as seen in previous work. How-ever, our results indicate that only slow coplanar encoun-ters (cases “Conf1” and “Conf2”), where the shock layerhad time to accrete substantial mass, can undergo frag-mentation. This differs from the results of Watkins et al.(1998b), where a large fraction of their objects formed inthe fragmentation of shock layers in both coplanar and non-coplanar cases. This discrepancy may be due to the use ofmore unstable initial disks, numerical problems related tolower resolution or the absence of vertical structure in thedisk models (Appendix A). However it should be noted thatcoplanar encounters are very rare if the disk orientations arerandom.The velocity constraint found in our simulations impliesthe number of BDs produced in encounters is a strong func-tion of the velocity dispersion in the young star cluster. How-ever, the numerical trend as a function of cluster density canbe complicated. On one hand, more massive and denser clus-ters usually have larger velocity dispersions and hence likelyto produce fewer sub-stellar objects per encounter. On theother hand, high number density will increase the encounterrate which could potentially increase the BD numbers.The encounter mechanism is complementary to othermechanisms for producing BDs. Proto-stellar disks in closeproximity are naturally produced when star cluster forma-tion is simulated from cloud collapse. Simulations of mas-sive clouds typically lack the resolution necessary to formlarge stable disks or to follow subsequent encounters in de-tail (Bate et al. 2003). Thus this mechanism is compatiblewith ideas such as competitive accretion and turbulent frag-mentation. The relative importance of these processes forforming proto-stars may influence the typical separation andorientation of proto-stellar disks and thus affect the startingpoint for encounter driven object formation. These statis-tics also depend the larger environment, such as the size,the mass and the amount of turbulent motion in the starforming cloud.A dynamical encounter origin for some BDs might ex-plain observed differences in binary statistics between starsand BDs in different cluster environments. Some observa-tions have indicated that the ratio of BDs to stars is higherin denser clusters. For example, Slesnick et al. (2004) foundin Trapezium, the number ratio ( R ss ) of sub-stellar objects(0 . < = M/M ⊙ < = 0 .
08) over stars (0 . < M/M ⊙ <
10) isabout 0.20, while Luhman et al. (2003) calculated for Tau-rus aggregates and IC 348 R ss is about 0.14 and 0.12 respec-tively, significantly lower than Trapezium. The results arenot conclusive as other surveys have indicated that R ss in Taurus is comparable with that in Trapezium (Levine et al.2006, and references therein). To make further predictionsabout how the BD ratio may be related to the cluster proper-ties, more comprehensive simulations and parameter studieswould be required, in which the encounter rate in differentcluster environment is estimated in detail and the encountervelocity parameter space is well sampled.Almost all the clumps formed in the isothermal simula-tions have flattened shapes. In the future evolution, we ex-pect these clumps to centrally condense to form a proto-BDwith a disk. Material should accrete onto the central objecton the viscous timescale. Observations indicate that the ac-cretion rate for low mass stars and sub-stellar objects canbe approximately fitted by ˙ M ∼ − M ⊙ yr − ( M/M ⊙ ) (Muzerolle et al. 2003). Hence for a relatively massive browndwarf with 0.05 M ⊙ the accretion rate is around 2 . × − M ⊙ yr − . With this rate the BD disk can havea mass as low as 10 − M ⊙ and still have lifetime com-parable to proto-stellar disks. To shed angular momen-tum, a proto-BD disk might also launch bipolar outflows orform sub-stellar companions. Since most BDs formed dur-ing encounter-induced disk fragmentation have large excessangular momenta, we expect outflows or companions to bebe common for young BD formed this way. Recent obser-vations have found forbidden emission lines in BD spec-tra, suggesting outflows (Natta et al. 2004). In particular,an outflow from a young BDs ρ Oph 102 has been spa-tially resolved (Whelan et al. 2005). Some of proto-BD disksin our simulations underwent secondary fragmentation andform binary BD systems (cf. Section 5.3). Due to resolu-tion and physical limitations, this result needs to be con-firmed with better simulations. However, the rapid rotationof these disk favours fragmentation (Matsumoto & Hanawa2003) which can form companions of a range of masses. Thisprovides a mechanism for forming the observed BD-BD bi-naries with the small separations peaking at ∼ R ss . Secondly, when the optical-thin approximation fails inhigh density regions, the effects of shock heating and ra-diative cooling become more important and should be fur-ther investigated. Including radiative transfer in the simu-lation is likely to reduce the number of objects formed perencounter, and provide for more realistic evolution of theproto-BD disks (as indicated in the optically thick approxi-mation, section 6). Alternately, with simplified models, theevolution of the clump themselves could be followed overlong timescales. c (cid:13) , 000–000 S. Shen et al.
APPENDIX A: RESOLUTIONREQUIREMENTS
When designing the simulation the resolution requirementwas estimated from the Bate & Burkert (1997) criterion,i.e., the Jeans mass should be resolved by at least twice thenumber of neighbour particles, N Neighbor . In our simulationswe used exactly 32 neigbours. In the initial condition, thedensity of our disk varies from 10 − g/cm to 10 − g/cm .Using our temperature profile, described in section 2, theJeans mass at the highest density for the initial disk wasestimated to be 0.0095 M ⊙ . Hence resolving the initial diskof 0.6 M ⊙ , which contains about 63 of these smallest Jeansmasses requires about 4000 particles for each disk. Duringthe encounter, however, as shown in Section 4, large inflowsof gas and unstable modes increase the density significantly.The density typically increases by three orders of magnitudelocally where fragmentation will then occur. Therefore to re-solve the smallest Jeans mass during the simulation signifi-cantly more particles were necessary. To correctly simulategas with densities up to 10 − g/cm about 127000 particleswere needed. We used 200000 particles for each disk, whichmet the above requirement and ensured that the fragmen-tation was physically modeled.The effects of resolution and the finite scale height werestudied with two-dimensional simulations and through an-alytical arguments by Nelson (2006). We performed twosets of comparison simulations to investigate these effectsin three dimensions. The first set consisted of three simu-lations of isolated disks with the same surface density andtemperature profile and a minimum Toomre Q of 1.3, wellabove the critical value of 0.8. The first disk was modeledby the single layer of 2000 particles. The second had 20000particles but still had no vertical structure. The third oneused 200000 particles and the vertical structure was mod-eled using the methods described in section 2. Figure A1depicts the results for each disk at a later time step.While no structure formed in the third, well-resolvedsimulation, strong gravitational instabilities and fragmenta-tion occurred in the first two within several outer orbitalperiods. Note that in the first case with 2000 particles theJeans mass was not resolved in the high density region of theinitial disk. For the second case although the Jeans mass wasfully resolved everywhere in the initial condition (not nec-essary during the simulation), the lack of vertical structurein the model was still destabilizing and artificial fragmenta-tion occurred. Although it is not clear from these two simu-lations which factor (resolving Jeans mass or modeling thefinite thickness) is more important to prevent unphysicalfragmentation, the test simulations show that it is essentialto consider both factors in the simulation.The second series consist of two simulations of retro-grade dynamical encounters of disks with Q = 2.1. The diskswere very stable in isolated runs. The dynamical parametersof the encounter was the same as case “Q3” in Table 1. Inthe first simulation the disks were modeled by single layerof 2000 particles, while in the second the disk was resolvedby 200000 particles and had vertical structures. The resultsis depicted in Figure A2.Similar to the isolated runs, without properly resolv-ing the Jeans mass and scale height, the encounter triggeredfragmentation at time step of 19700 years in the shock layer between two coplanar disks (Panel 2, Figure A2). In the sec-ond simulation, although the shock layer was still produced(Panel 4), no fragmentation occurred in the layer. The en-counter was also modeled using 400000 particles with thesame result.In summary, regardless of whether the disks are pas-sively evolving or undergo a physical impact, numerical ef-fects may be significant if the Jeans mass is not resolved orif no vertical structure is presented in the disk model. Themajor numerical effect is then enhanced instability and theformation of artificial clumps. Our resolution of 200000 par-ticles was chosen based on a conservative estimate of resolu-tion requirements, e.g., it assumed the density of the wholedisk could increase up to 10 − g/cm . In fact, most outerregion of the disk had density below 10 − g/cm when frag-mentation initial occurred. With this resolution it was en-sured that the initial disk fragmentation (if any) was phys-ical. We also performed simulations at higher resolution(400000 and 800000 particles) for some specific cases andthe results converged. REFERENCES
Andre, P., Ward-Thompson, D., & Barsony, M. 2000, Pro-tostars and Planets IV, 59Barrado y Navascu´es, D., Stauffer, J. R., Bouvier, J.,Jayawardhana, R., & Cuillandre, J.-C. 2004, ApJ, 610,1064Bate, M. R., & Burkert, A. 1997, MNRAS, 288, 1060Bate, M. R., Bonnell, I. A., & Bromm, V. 2002, MNRAS,332, L65Bate, M. R., Bonnell, I. A., & Bromm, V. 2003, MNRAS,339, 577Boffin, H. M. J., Watkins, S. J., Bhattal, A. S., Francis, N.,& Whitworth, A. P. 1998, MNRAS, 300, 1189Boss, A. P. 2001,ApJ, 563, 367Brice˜no, C., Luhman, K. L., Hartmann, L., Stauffer, J. R.,& Kirkpatrick, J. D. 2002, ApJ, 580, 317Burrows, A., et al. 1998, Brown Dwarfs and ExtrasolarPlanets, 134, 354Caballero, J., Bjar, V., Rebolo, R., Eisloeffel, J., ZapateroOsorio, M., Mundt, R., Barrado Y Navascus, D., Bihain,G., Bailer-Jones, C., Forveille, T., Martn, E. 2007, A&A,47[0, 903Chabrier, G. 2002, ApJ, 567, 304Chauvin, G., Lagrange, A.-M., Dumas, C., Zuckerman, B.,Mouillet, D., Song, I., Beuzit, J.-L., & Lowrance, P. 2005,A&A, 438, L25Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368D’Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553,321Elmegreen, B. G. 1999, ApJ, 522, 915Gammie, C. F. 2001, ApJ, 553, 174Gauvin, L. S., & Strom, K. M. 1992, ApJ, 385, 217Gizis, J. E., Kirkpatrick, J. D., Burgasser, A., Reid, I. N.,Monet, D. G., Liebert, J., & Wilson, J. C. 2001, ApJL,551, L163Gorlova, N. I., Meyer, M. R., Rieke, G. H., & Liebert, J.2003, ApJ, 593, 1074Guieu, S., Dougados, C., Monin, J.-L., Magnier, E., &Mart´ın, E. L. 2006, A&A, 446, 485 c (cid:13) , 000–000 ncounter-induced brown dwarf formation −13 −6−13 −6−13 −7t = 9,000 years t = 6, 000 years t = 24, 000 years1 2 3 Figure A1.
Isolated disk simulations with different resolutions and disk models. All three disks have Q = 1.3 initially. Time steps foreach snapshot are given in the upper left corners. Boxed regions indicate the fragments (if any). The grey scale under each panel givesthe density ranges in units of M ⊙ /AU . −13 −6 −10 −6 −13 −6 −10 −6t = 19.700 yearst = 19,700 years t = 21,000 years t = 21,000 years1 2 43 Figure A2.
Disk encounter simulations with different resolutions and disk models. The disks have Q = 2.1 initially. Time steps for eachsnapshot are given in the upper left corners. The first two panels (panel 1-2) depicted the low resolution run and panel 3-4 are from thehigher resolution one. Boxed regions indicate the fragments (if any). The grey scale under each panel gives the density ranges in units of M ⊙ /AU . Hubber, D. A., & Whitworth, A. P. 2005, A&A, 437, 113Jayawardhana, R., Mohanty, S., & Basri, G. 2003, ApJ,592, 282Johnson, B. M., & Gammie, C. F. 2003, ApJ, 597, 131Kim, W.-T., Ostriker, E. C., & Stone, J. M. 2002, ApJ,581, 1080Kraus, A. L., White, R. J., & Hillenbrand, L. A. 2006, ApJ,649, 306Kroupa, P. 2001, MNRAS, 322, 231Krumholz, M. R., McKee, C. F., & Klein, R. I. 2005, Na-ture, 438, 332Levine, J. L., Steinhauer, A., Elston, R. J., & Lada, E. A.2006, ApJ, 646, 1215Lin, D. N. C., Laughlin, G., Bodenheimer, P., & Rozyczka,M. 1998, Science, 281, 2025Lodato, G., Meru, F., Clarke, C. J., & Rice, W. K. M. 2007,MNRAS, 374, 590Luhman, K. L. 2000, ApJ, 544, 1044Luhman, K. L., Stauffer, J. R., Muench, A. A., Rieke, G. H., Lada, E. A., Bouvier, J., & Lada, C. J. 2003, ApJ,593, 1093Mart´ın, E. L., Dougados, C., Magnier, E., M´enard, F.,Magazz`u, A., Cuillandre, J.-C., & Delfosse, X. 2001, ApJ,561, L195Matsumoto, T., & Hanawa, T. 2003, ApJ, 595, 913Marcy, G. W., & Butler, R. P. 2000, PASP, 112, 137Mayer, L., Quinn, T., Wadsley, J., & Stadel, J. 2004, ApJ,609, 1045Mejia, A. C. 2004, Ph.D. Thesis,Mohanty, S. 2005, Protostars and Planets V, 8532Moraux, E., Bouvier, J., Stauffer, J. R., & Cuillandre, J.-C.2003, A&A, 400, 891Muench, A. A., Alves, J., Lada, C. J., & Lada, E. A. 2001,ApJL, 558, L51Mundy, L. G., Looney, L. W., & Welch, W. J. 2000, Pro-tostars and Planets IV, 355Muzerolle, J., Hillenbrand, L., Calvet, N., Brice˜no, C., &Hartmann, L. 2003, ApJ, 592, 266 c (cid:13) , 000–000 S. Shen et al.
Natta, A., Testi, L., Muzerolle, J., Randich, S., Comer´on,F., & Persi, P. 2004, A&A, 424, 603Nelson, A. F., 2006, MNRAS, 373, 1039Padoan, P., & Nordlund, ˚A. 2004, ApJ, 617, 559Pickett, B. K., Mej´ıa, A. C., Durisen, R. H., Cassen, P. M.,Berry, D. K., & Link, R. P. 2003, ApJ, 590, 1060Price, D. J. & Bate, M. R. 2009, NMNRAS, accepted
Rafikov, R. R. 2007, ApJ, 662, 642Reipurth, B., & Clarke, C. 2001, AJ, 122, 432Rice, W. K. M., Lodato, G., & Armitage, P. J. 2005, MN-RAS, 364, L56Simon, M., Bender, C., & Prato, L. 2006, ApJ, 644, 1183Slesnick, C. L., Hillenbrand, L. A., & Carpenter, J. M. 2004,ApJ, 610, 1045Stamatellos, D. & Whitworth, A. P. 2009, MNRAS, 392,413Thies, I., Kroupa, P., & Theis, C. 2005, MNRAS, 364, 961Toomre, A. 1964, ApJ, 139, 1217Toomre, A., & Toomre, J. 1972, ApJ, 178, 623Truelove, J. K., Klein, R. I., McKee, C. F., Holliman, J. H.,II, Howell, L. H., & Greenough, J. A. 1997, ApJL, 489,L179Wadsley, J. W., Stadel, J., & Quinn, T. 2004, New Astron-omy, 9, 137Watkins, S. J., Bhattal, A. S., Boffin, H. M. J., Francis, N.,& Whitworth, A. P. 1998, MNRAS, 300, 1205Watkins, S. J., Bhattal, A. S., Boffin, H. M. J., Francis, N.,& Whitworth, A. P. 1998, MNRAS, 300, 1214Whelan, E. T., Ray, T. P., Bacciotti, F., Natta, A., Testi,L., & Randich, S. 2005, Nature, 435, 652White, R. J., Greene, T. P., Doppmann, G. W., Covey,K. R., & Hillenbrand, L. A. 2007, Protostars and PlanetsV, 117Whitworth, A., Bate, M. R., Nordlund, ˚A., Reipurth, B.,& Zinnecker, H. 2007, Protostars and Planets V, 459Yorke, H. W., & Bodenheimer, P. 1999, ApJ, 525, 330Vorobyov, E. & Basu, S., 2009, MNRAS, 393, 822 c (cid:13)000