Spiral Arms in Disks: Planets or Gravitational Instability?

Spiral arm structures seen in scattered light observations of protoplanetary disks can potentially serve as signposts of planetary companions. They can also lend unique insights into disk masses, which are critical in setting the mass budget for planet formation but are difficult to determine directly. A surprisingly high fraction of disks that have been well-studied in scattered light have spiral arms of some kind (8/29), as do a high fraction (6/11) of well-studied Herbig intermediate mass stars (i.e., Herbig stars $>1.5M_\odot$). Here we explore the origin of spiral arms in Herbig systems by studying their occurrence rates, disk properties, and stellar accretion rates. We find that two-arm spirals are more common in disks surrounding Herbig intermediate mass stars than are directly imaged giant planet companions to mature A and B stars. If two-arm spirals are produced by such giant planets, this discrepancy suggests that giant planets are much fainter than predicted by hot start models. In addition, the high stellar accretion rates of Herbig stars, if sustained over a reasonable fraction of their lifetimes, suggest that disk masses are much larger than inferred from their submillimeter continuum emission.As a result, gravitational instability is a possible explanation for multi-arm spirals. Future observations can lend insights into the issues raised here.


Introduction
As the birthplaces of planets, protoplanetary disks surrounding young stars lend unique insights into the circumstances under which planets form. The initial masses of disks set an upper limit on the mass budget for planet formation, a basic constraint on theories of how planets form. From the detailed study of disk morphologies, we may also glean evidence of the formation of planets themselves, e.g., from the detection of gaps, rings, and spiral arms created by planets. As a result, numerous studies have examined the protoplanetary disks around T Tauri stars and Herbig stars for clues to their structure and planet formation status.
Despite their fundamental importance, disk masses have been proven difficult to measure. Although disk masses are commonly inferred from submillimeter continuum fluxes, which probe the dust content of disks, several arguments strongly suggest that much of the solids in T Tauri disks have grown beyond the sizes probed by submillimeter continua. As one simple example, submillimeter continuum disk masses are too low to even account for the solids that are known to be locked up in exoplanet populations and debris disks, indicating that T Tauri disks harbor more massive reservoirs of solids than those probed by continuum measurements (Najita & Kenyon 2014), and disk masses are underestimated as a result.
More graphically, high resolution submillimeter continuum images of disks made with ALMA commonly show stunning rings and gaps (ALMA Partnership et al. 2015;Andrews et al. 2016;Cieza et al. 2016;Isella et al. 2016;Cieza et al. 2017;Loomis et al. 2017;Fedele et al. 2018;Huang et al. 2018;Dipierro et al. 2018), structures that signal that disk solids have grown beyond centimeter sizes and possibly into planetary mass objects (e.g., Dong et al. 2015c;Dipierro et al. 2015;Dong et al. 2017a;Bae et al. 2017). These results exclaim that disks are highly evolved and call into question disk mass estimates made under the usual assumption of simple, unevolved, optically thin disks. Given the difficulty of measuring disk masses through proxies such as dust continuum emission, other approaches are welcome.
Alongside these developments, high contrast images of disks made in scattered light also reveal remarkable features. Near-infrared scattered light imaging observations probe the spatial distribution of small grains, approximately micron-sized or smaller, that are wellcoupled to the gas (Zhu et al. 2012;Pinilla et al. 2012). These studies have revealed spiral arms in a growing number of disks (Figure 1), both prominent near-symmetric two-arm spirals (e.g., Grady et al. 2013;Benisty et al. 2015) and more flocculent multi-arm spirals on smaller scales and at lower contrast (e.g., Fukagawa et al. 2004;Hashimoto et al. 2011). Spiral arms can be produced by either gravitational instability (e.g., Rice et al. 2003;Lodato & Rice 2004;Stamatellos & Whitworth 2008;Kratter et al. 2010;Kratter & Lodato 2016) or the presence of massive companions (gas giant to stellar masses; e.g., Kley & Nelson 2012;Zhu et al. 2015;Muñoz & Lai 2016).
The presence or absence of spiral arms in disks, if induced by gravitational instability, can place a dynamical constraint on disk masses. Gravitational instability occurs when the destabilizing self-gravity of a disk dominates over the restoring forces of gas pressure and differential rotation (Toomre 1964;Goldreich & Lynden-Bell 1965), a requirement that roughly translates into M disk /M 0.1 under typical conditions (Kratter & Lodato 2016). Spiral arms induced by gravitational instability are expected to appear nearly symmetric in scattered light imaging, with the number of arms roughly scaling as n ≈ M /M disk (Lodato & Rice 2004Cossins et al. 2009;Dong et al. 2015a).
Planetary companions can also drive spiral arms (Goldreich & Tremaine 1979;Ogilvie & Lubow 2002). A massive companion (above a few × 10 −3 M ) generates a set of nearsymmetric two-arm spirals interior to its orbit, similar to the morphology of observed twoarm spirals (Dong et al. 2015b, see also Fig. 2a,b), with the companion located at the tip of the primary arm (e.g., Zhu et al. 2015;Fung & Dong 2015;Bae et al. 2016). Lower mass companions (∼ 10 −3 M ) may generate additional arms interior to their orbits, weaker than the primary and secondary arms in scattered light (e.g., Juhász et al. 2015;Fung & Dong 2015;Lee 2016;Bae et al. 2017;Bae & Zhu 2018, see also Fig. 2d,e). If spiral arms are driven by planets, arm morphologies can provide clues to the location and masses of the planets.
Here we explore these hypotheses for the origin of the spiral arms detected in scattered light by comparing (i) the demographics of arm-bearing disks (disk masses and stellar accretion rates) with those of other young stars and (ii) the incidence rate of arms compared to that of resolved planetary companions to mature stars. Although we begin by considering the detection rates of spiral arms in disks surrounding both T Tauri stars and Herbig stars, we end up focusing more closely on the properties of the higher mass Herbig stars, where the census of scattered light structures is more complete.
In §2 we characterize the occurrence rate of spiral arm structures in scattered light imaging of disks. To explore whether spiral arms could be produced by orbiting companions, we compare in §3 the occurrence rates of spiral arms in disks to that of giant planetary companions to mature stars. In §4 we explore the alternative possibility that spiral arms are driven by gravitational instaibility by examining the masses of disks that show spiral arm structures, as probed by submillimeter continuum emission and stellar accretion rates. In §5 we examine the planet and gravitational instability hypotheses in greater detail before suggesting scenarios that are potentially consistent with our results. We summarize our results in §6, where we also describe future studies that can lend insights into the issues raised here.

Sample Selection and Statistics
In near-infrared scattered light imaging of disks, spiral arms have been found in both embedded Class 0/I objects surrounded by infalling envelopes (Canovas et al. 2015;Liu et al. 2016b) and revealed Class II objects (i.e., Herbig and T Tauri stars; Figure 1). While gravitational instability may drive the arms in Class 0/I disks (e.g., Vorobyov & Basu 2005Tobin et al. 2016;Pérez et al. 2016;Dong et al. 2016b;Tomida et al. 2017;Meru et al. 2017), which tend to be massive and are fed by infall from the envelope, the origin of the arms detected in Class II disks is unclear. Figure 2 shows how in principle spiral arms could serve as signposts of either an orbiting planet or a gravitationally unstable disk. The left column shows observed images of two representative disks: MWC 758, with its prominent near-symmetric two-arm spirals (Fig. 2a); and HD 142527, with multiple small-scale, lower contrast arms (Fig. 2d). The middle and right columns show synthetic observations from the literature, produced in combined hydrodynamic and radiative transfer simulations, which illustrate how morphologically similar spiral structures can occur in disks with orbiting planets (middle column) or when a disk is gravitational unstable (right column). As discussed in §1, both a multi-Jupiter mass giant companion and a gravitationally unstable disk with M disk ∼ 0.5M can produce a set of prominent two-arm spirals (Fig. 2b,c), while lower mass planets (M p ∼ M J ) and marginally gravitationally unstable disks with M disk ∼ 0.1M can produce multiple small-scale arms at low contrast (Fig. 2e, f).
To attempt to distinguish between these two scenarios for the origin of spiral structure in disks, we first characterize, in this section, the occurrence rate of spiral arms in scattered light imaging of disks. §2.1 describes the data available and our process of selecting a sample of disks well studied in scattered light that are suitable for our study. The occurrence rate of arms within this sample is described in §2.2. Because the sample of well-studied disks is small in size, in §2.3 we examine the occurrence rate of arms within in a volume-limited sample of Herbig stars, which provides a lower limit on the true occurrence rate.

Disk Imaging Sample
We collected from the literature a sample of protoplanetary disks that have been spatially resolved in scattered light at optical to NIR wavelengths (Appendix Table A .1.1) and that meet four main criteria. First, the disk is gas-rich (i.e., debris disks are excluded). Second, the disk was observed with an inner working angle < 1 , small enough to reveal the inner ∼100 AU planet forming region (most disks imaged to date are 100-200 pc away).
Third, the disk surface was revealed with sufficient signal-to-noise to identify possible structures. Finally, the disk was observed with an angular resolution 0. 1 in order to detect and resolve spiral arms.
The third condition largely restricts the sample to Class II objects without a substantial envelope, as the disk surface of more embedded sources is often not accessible in scattered light. It also excludes disks that are completely flat or shadowed (Dullemond & Dominik 2004). The last criterion restricts the sample to those observed at L-band (3.8 µm) or shorter wavelengths, because λ/D 0. 1 requires λ 3.8 µm for observations made with D ≈ 8 m diameter telescopes. Systems that have been directly imaged, both with and without a coronagraphic mask, are included, whereas interferometric observations are not. The resulting sample of 49 objects (Table A.1.1) derives mainly from a few major exoplanet and disk surveys, such as the Subaru-based SEEDS (Tamura 2009;Uyama et al. 2017) and the VLT-based SHINE surveys, supplemented by studies of smaller samples that used a variety of ground-based instruments and HST. Roughly half of the sources are transitional disks (i.e., disks whose inner region is optically thin in the continuum; Espaillat et al. 2014).
We further removed sources viewed at high inclination ( 70 • ) whose scattered light structures may be difficult to assess (AK Sco, T Cha, AA Tau, RY Lup, HH 30, HK Tau, HV Tau C, LkHα 263C, and PDS 144N). It is possible to identify spiral arms in disks at such high inclinations (e.g., RY Lup; Langlois et al. 2018); however it is usually a challenge because of the resulting geometric effect (Dong et al. 2016a). We also excluded any remaining Class 0/I sources and sources undergoing accretion outbursts (FU Ori, HL Tau, R Mon, V 1057 Cyg, V1735 Cyg, and Z CMa) in order to probe a restricted evolutionary state and to avoid scattered light contamination from an infalling envelope.
Sources with known stellar companions within 5 (HD 100453, GG Tau, HD 150193, and T Tau) were also excluded to avoid systems in which a stellar companion induces disk structure. Whereas a massive planet exterior to the arms can generate the observed two-arm spirals, a stellar companion closer than ∼ 5 could produce similar spiral structure in a disk with a typical size of ∼1−2 ; (e.g., as in the case of the HD 100453 disk; Wagner et al. 2015;Dong et al. 2016c;Wagner et al. 2018). The companion to HD 142527, which is located well within the gap in the disk, is commonly interpreted as having too low a mass (0.13M ; Lacour et al. 2016) and eccentricity to generate the observed spiral arms, so this object is included in our sample (although see Price et al. 2018). Lastly, we removed HD 141569 from the sample, because its low dust mass (on the order of 1 M ⊕ or less; Flaherty et al. 2016;White et al. 2016) makes the system prone to radiation pressure and photoelectric instability, which can also produce spiral arm structures (e.g., ?).
Finally, we arrive at a sample of 29 "well-studied NIR disks", whose stellar and disk -6properties are listed in Table 1. The tabulated M values were obtained from T eff and L using the Siess et al. (2000) pre-main-sequence evolutionary tracks. Notes on individual sources are provided in Appendix §D. The final sample contains 11 Herbig sources (the first 11 rows in the table) and 18 T Tauri stars. Of the 29 sources, 8 have spiral arms: the MWC 758, SAO 206462, LkHα 330, and DZ Cha disks show prominent two-arm spirals ( Fig. 1, top), the AB Aur, HD 142527, and HD 100546 disks show multiple small arms, and the V1247 Ori disk has one arm ( Fig. 1, bottom).

Spiral Arm Statistics Among Well-studied NIR Disks
From the results shown in Table 1, we find that the incidence rate of spiral arms in disks is surprisingly high among young well-studied A and B stars (3/7), well-studied FGK stars (5/22), and well-studied disks overall (8/29; Table 3).
We can make a similar accounting of spiral fraction by stellar mass rather than by spectral type. Anticipating our comparison of the spiral arm incidence rate to that of massive companions to mature AB stars ( §3), we define a Herbig intermediate mass star (IMS) sample over a broader range of spectral type. A main sequence F0 star has a mass of ∼ 1.5M (e.g., Pecaut & Mamajek 2013). So we consider as Herbig IMS all Herbig systems with spectral types between F6 and B9, and stellar masses > 1.5M . Table 3 shows that the spiral fraction is 6/11 among well-studied Herbig IMS and 2/18 among the low mass stars (LMS) that make up the rest of the sample. The incidence rate of two-arm spirals is also high: 2/11 among the well-studied Herbig IMS and 2/18 among the LMS. Although the samples are small, the data can nevertheless lend insight into whether there is a significant difference between the occurrence rates of spiral structure in LMS and Herbig IMS disks. For example, using Fisher's exact test (a contingency table analysis), we find that there is a 3% probability (two-tailed p-value) that the spiral occurrence rate among Herbig IMS is the same as that among low mass stars. Thus the apparent low occurrence rate of spirals among low mass stars compared to Herbig IMS is statistically significant (see also Avenhaus et al. 2018;Isella & Turner 2018).
There are a couple of biases associated with these results. Firstly, scattered light imaging studies often target sources with interesting known structures (e.g., transitional disks) and sources that appear to show interesting structures at modest signal-to-noise receive more extensive follow up. Secondly, scattered light imaging studies are only sensitive to disk structures that can be seen in scattered light. If a disk's surface is not sufficiently flared, its spiral structure cannot be detected with this technique.
-7 -Despite these difficulties, we can obtain a lower limit on the true spiral fraction of disks by examining the fraction of disks with known spiral arms within a volume-limited sample of comparable young stellar objects. We therefore focus our study on the Herbig IMS population, because a fair fraction of such sources within the local volume have already been surveyed in scattered light. In contrast, only a small fraction of the (lower mass) T Tauri population has been studied, and as a result, a corresponding lower limit on the true spiral fraction for LMS would not be a significant constraint.

A Volume-limited Sample of Herbig Stars
Excluding V1247 Ori, which is very distant (320 pc), the remaining Herbig IMS in Table 1 are all within 200 pc. We therefore focus on a sample of Herbig IMS within this volume, i.e., Herbig sources within 200 pc that have spectral type F6 or earlier, M > 1.5M , and no known stellar companion within 5 . Table 2 compiles a list of sources from the literature meeting these criteria The et al. 1994;Malfait et al. 1998;Vieira et al. 2003;Erickson et al. 2011;Chen et al. 2012). Generally, these stars are identified by their association with star forming regions or a reflection nebula and the presence of hydrogen emission lines (particularly Hα). Because the sample could be useful for studies other than the one carried out here, the list also includes, for completeness, sources that we ignore in our study (indicated in Column 4 of Table 2): sources with stellar companions 0. 3-5 away (HR 811, V892 Tau, PDS 178, CQ Tau, HD 100453, HR 5999, MWC 863, TY CrA, T CrA), sources with M < 1.5M (AK Sco), sources with low, potentially optically thin dust (HD 141569), and sources with decretion disks (51 Oph).
Excluding these sources, there are 24 remaining sources, which we carry forward as our "volume-limited Herbig IMS" sample. Ten of these have been well studied in NIR imaging (i.e., they are in Table 1). Note that with the above prescription for the Herbig IMS sample, we exclude younger sources such as LkHα 330 (Table 1) which also has a stellar mass > 1.5M but a later spectral type (G3; Figure 3), These sources are excluded in part because earlier spectral type sources (F5-B9) have been studied more completely in scattered light; younger, more embedded disks are more difficult to study in scattered light due to contamination from the envelope. They are also more numerous and much less completely studied in scattered light. The masses of the later spectral type sources are more uncertain as well.
Within the volume-limited Herbig IMS sample ( Table 2), 5 of the 10 well-studied sources display spiral structure (AB Aur, MWC 758, HD 100546, SAO 206462, HD 142527). The -8true spiral fraction depends on whether the remaining 14 sources that have not been well studied (i.e., are not included in Table 1) have arms or not. Thus, the minimum spiral fraction is 5/24, and the minimum two-arm spiral fraction is 2/24 (Table 4). We conclude that even after correcting for systems that have not been studied well in scattered light, the spiral fraction of Herbig IMS is high. Detailed studies of the additional Herbig IMS within 200 pc would be extremely valuable in firming up the spiral fraction.

Comparison with Exoplanet Demographics
The two-arm spirals observed in scattered light can be produced by 5-13 M J planets located 30-300 AU from the star Dong 2015 andFung 2017). To explore whether this scenario is tenable, we compare the incidence rate of arms with the incidence rate of planetary companions to mature A and B stars in the corresponding mass and separation range. Planetary companions to mature stars are an advantageous population to study because they are better characterized than are companions to Herbig stars and T Tauri stars (Bowler 2016).
Compiling the results of direct imaging searches for planetary companions in this mass and separation range, Bowler (2016) reported that only 3 mature (single) AB stars have such a companion out of a total sample of 110 stars with spectral types B2 to A9. This result assumes that planetary companions have the properties predicted by "hot start" models (see §5). In comparison, the two-arm spiral fraction among well-studied Herbig IMS is larger (2/11; §2.2). Although the number of detections is small in both cases, there is only a 6% probability of obtaining such discrepant results if the the two rates are the same, according to Fisher's exact test.
Further observations are needed to explore this possible difference. As noted in §2.3), the two-arm spiral fraction in the volume-limited Herbig IMS sample (Table 4) is at least 2/24. If future observations show that there are no other two-arm spirals within 200 pc than the two already known, the probability that the two rates are the same would be 21% and the difference in the observed rates would be of little statistical significance. However, if two additional two-arm spirals are found in the sample, there would only be a 2% probability that the 4/24 spiral fraction and the Bowler (2016) fraction are the same, and the difference in the rates would be quite significant.
The occurrence rate of companions to mature FGK stars is lower than that for mature AB stars (Bowler 2016). To compare with the FGK companion rate, we can select sources from Table 1 that are "future FGK stars", i.e., sources with stellar masses between 0.5 M and 1.5 M (TW Hya, PDS 70, J1604-2130, RX J1615.3-3255, DoAr 28, V4046 Sgr, GM Aur, LkCa 15, PDS 66, SR 21, IM Lup, DZ Cha). The incidence rate of two-arm spirals among the well-studied "future FGK stars" is 1/12 compared with the 0/155 incidence rate of 5-13 M J companions to mature FGK stars in the 30-300 AU separation range (Bowler 2016;Vigan et al. 2017). There is a 7% probability that the two rates are the same, similar to the situation for the higher mass stars. Further observations are needed to quantify the extent to which these rates differ.
With the apparent paucity of massive planetary companions to mature AB stars appearing to question the idea that such planetary companions drive two-arm spirals in Herbig disks, in the next section we explore the alternative possibility that spiral arms are driven by gravitational instability. We return in §5 to the planetary companion hypothesis and discuss ways in which it may be consistent with observations.

Disk Masses and Gravitational Instability
To help interpret the origin of the spiral arms detected in scattered light imaging, we have also compiled disk masses and stellar accretion rates for the volume-limited Herbig IMS sample ( Table 2). The latter can provide a complementary estimate of disk (gas) mass assuming a typical accretion history for the disk (Hartmann et al. 1998). Tables 1 and 2 tabulate the disk dust masses M dust calculated from submillimeter continuum fluxes (Appendix §B). Stellar accretion ratesṀ are determined from the veiling of the Balmer discontinuity or hydrogen emission line fluxes; further details are provided in Appendix §C.
Turning first to the dust masses, although they may not be a reliable quantitative measure of disk mass ( §1), they might serve as a relative measure: sources with higher disk gas masses may also have higher dust masses. Assuming M disk = 100M dust , Figure 4 shows the ratio of M disk /M for all sources in the volume-limited Herbig IMS sample that have estimated dust masses. The ratio M disk /M ranges from ∼ 10 −3 − 10 −1 . Apart from HD 142527, which has a very large dust mass, the sources with two arms (blue bars) or multiple arms (green bars) are distributed throughout this range and interspersed with other sources, both those known to have no arms (black bars) and those that have not yet been studied in scattered light in detail (gray bars).
There is no obvious difference in the M disk /M ratio of disks with and without arms. The average ratio is 0.008 for the entire volume-limited Herbig IMS sample. Sources with arms and those known to have no arms both have the same average ratio, 0.02. In other words, scattered light surveys have preferentially studied sources with higher submillimeter -10continuum fluxes.
One indication that these values are flawed as quantitative measures of disk mass comes from estimating the remaining disk lifetime assuming these disk masses, τ life = M disk /Ṁ . Figure 5 shows τ life for sources in the volume-limited Herbig IMS sample for which both measurements are available. If we assume that the remaining disk lifetime is 2 Myr, comparable to the average age of disk-bearing young stars (Hernández et al. 2008;Richert et al. 2018) and a conservative estimate based on the HR diagram (Figure 3), half of the sources have lifetimes 0.2 Myr, a mere 10% of the average age. The short lifetime is a consequence of the high accretion rates. If these values are correct, half of the Herbigs within 200 pc are in the last 10% of their mass-building phase of life, which seems implausible.
We could instead assume that we are viewing Herbigs on average at middle age, a current age of at least t 0 ∼ 2 Myr. Parametrizing the decline ofṀ with time asṀ (t) = M (t 0 )(t/t 0 ) −η (e.g., Hartmann et al. 1998), where t 0 is the current age of the star, the total mass that the star accretes at later times -Aguilar et al. 2010, Hartmann et al. 1998 13, about ten times higher than the disk masses estimated from submillimeter continua. If the stars are older, the implied disk masses are higher. Figure 6 shows M disk /M where M disk is estimated as above. The sources with two (blue bars) or multiple arms (green bars) are concentrated in the upper half of the sample, reaching values as large as 0.6 (AB Aur), with about one-half to one-third of the sample above the nominal gravitational instability limit of M disk /M = 0.1. The clustering of "arm" sources toward higher M disk /M would be consistent with the arms being generated through or enhanced by gravitational instability. It would extremely interesting to measure the scattered light morphology of the unstudied sources (gray bars). If the lower mass disks also show spiral arms, it is unlikely that gravitational instability plays the major role in generating arms. Figure 7 compares the above submillimeter continuum-based and stellar accretion-based disk masses for the volume-limited Herbig IMS sample as well as the T Tauri stars studied by Andrews & Williams (2007). Both groups of sources have, on average, 3-10 times larger masses estimated from stellar accretion rates than from dust continuum. The discrepancy between the disk masses estimated by these two methods has been previously noted for T Tauri stars by several authors (e.g., Andrews & Williams 2007;Rosotti et al. 2017). The Herbig IMS studied here show a similar trend.

Discussion
Among disks that have been well-studied in scattered light, spiral arm structures are surprisingly common ( §2.2), although their origin is unclear. While giant planets and gravitational instability have both been put forward as explanations for the observed spiral structure, neither explanation is readily confirmed in our study. We find that despite the small sample sizes studied to date, we can infer that two-arm spirals occur more frequently among Herbig IMS disks than do massive giant planet companions to mature AB stars, challenging the planet hypothesis ( §3). The viability of gravitational instability as an explanation for spiral structure is uncertain, because submillimeter continuum-based disk masses are low. However, as we showed, disk masses may be substantially larger if stellar accretion rates are a guide ( §4). In this section, we examine the proposed hypotheses in greater detail before presenting scenarios that are potentially consistent with our results.
Shadows and Stellar Flybys. Spiral arms detected in scattered light have also been modeled as arising from rotating shadows cast by a warped inner disk (e.g., Montesinos et al. 2016; ?) and stellar flybys (e.g., Pfalzner 2003). These alternative explanations require specific conditions that are unlikely to be satisfied in the majority of cases. The moving shadow scenario is best at explaining grand design two-arm spirals; yet in the only disk with both a clear m=2 shadow pattern and a set of two-arm spirals (HD 100453), the disk's arms are almost certainly induced by the visible stellar companion (Dong et al. 2016c;Wagner et al. 2018). To catch the spirals induced by a stellar flyby, one needs to image the system within a few dynamical times (i.e., < 10 4 years at typical distances) after the closest encounter before the spirals winded up due to differential rotation (Pfalzner 2003). The short timescales make it unlikely to catch the disk in this situation. Because these two options are very unlikely to be the general mechanisms responsible for observed spiral arms in scattered light, we focus more closely on the gravitational instability and planetary companion hypotheses that have been considered more intently in the recent literature.
Gravitational Instability. Exciting two-arm spirals in a gravitationally unstable disk requires M disk /M ∼ 0.5, a situation that may occur in Herbig disks, especially if they have recently decoupled from their infalling envelopes. In such systems, spiral shocks transport angular momentum and drive rapid accretion, withṀ 10 −5 M yr −1 (Kratter et al. 2008(Kratter et al. , 2010Dong et al. 2015a). Consequently, the disk mass decreases quickly and the two-arm phase is brief (on the order of 10 4 yr). As the disk mass drops, the disk accretion rate also declines, and the disk structure shifts to a multi-arm morphology (e.g., Rice & Armitage 2009, Hall et al, in prep.). Eventually the disk is stabilized against gravitational instability and the disk settles down to an even lower accretion rate.
This scenario appears inconsistent with the demographic properties of disks with two--12arm spirals. In the above picture, two-arm spirals should occur in systems with the largesṫ M . Instead, observed two-arm spiral disks have middle-of-the-roadṀ values among the Herbig IMS within 200 pc ( Figure 5). Two-arm spirals are also expected to be rare, because this phase lasts for on the order of 1% of the 2 Myr average age of disk-bearing stars once envelope infall have ceased. In contrast, two-arm spirals are found in 2/11 (or ∼18%) of the well-studied Herbig IMS and ≥ 2/24 (or 8%) of the volume-limited Herbig IMS sample ( Table 4). Disks with two-arm spirals are also expected to be evolutionarily young, and to have just emerged from their infalling envelopes. In the HR diagram, disks with two-arm spirals are of preferentially later spectral type than other well-studied Herbig stars, but they are not preferentially younger (Fig. 3).
Although it seems unlikely to explain the two-arm spirals, gravitational instability is a possible explanation for multi-arm spirals. Disks with M disk /M ∼ 0.1 have multiple arms and are marginally unstable, with disk accretion rates ∼ 10 −8 -10 −6 M yr −1 in the outer disk (e.g., Dong et al. 2015a;Hall et al. 2016). They can remain in this state for a period τ ∼ M disk /Ṁ , or 10 5 -10 7 yr, comparable to the 2 Myr average age of disk-bearing stars. Additional disk accretion mechanisms, such as the magnetorotational instability (Balbus & Hawley 1992;Gammie 2001) or magnetized disk winds (e.g., Bai et al. 2016;Bai 2017), will become more dominant as the gravitational instability-driven accretion rate declines.
The relatively long lifetime of the multi-arm phase of gravitational instability is possibly consistent with the 3/11 (or 27%) incidence rate of multi-arm spirals among well-studied Herbig IMS and ≥ 3/24 (or ≥ 13%) incidence rate within the volume-limited Herbig IMS sample (Table 4). All three of the multi-arm disks have large disk masses based on their stellar accretion rates (M disk /M 0.1; Figure 6) as well as stellar accretion rates in the above range ∼10 −8 -10 −6 M yr −1 ( Figure 5). Note that the similarity between observed stellar accretion rates and the predicted outer disk accretion rates suggest that gravitational instability-driven accretion in the outer disk could potentially be sustained down to small disk radii and onto the star; further theoretical work is needed to understand the details of the accretion process. Numerical simulations that quantify the time that a disk stays in each n-arm state in a realistic environment will be extremely useful to explore this possibility.
Planetary Companions. By studying sources that have been well studied in scattered light, we have focused on older, revealed Herbig stars that are not obscured by an infalling envelope. Thus our sample is biased against gravitational instability, which is more likely to occur in (younger) massive disks, and possibly toward giant planets, which may take a few Myr to form. It is therefore surprising to find that the properties of the sample seem inconsistent with the reported properties of giant planets around mature A stars: two-arm spirals appear to be more common than the giant planets massive enough to drive them.
-13 -One possible resolution to this discrepancy is that giant planets are fainter than we think, i.e., current direct imaging searches for planetary companions to mature A stars may underestimate the planetary masses to which they are sensitive. In converting companion fluxes (and detection upper limits) to masses, one can choose between "hot start" (e.g., COND and DUSTY models, Baraffe et al. 2003;Chabrier et al. 2000) and "cold start" evolutionary models (e.g., Fortney et al. 2008;Spiegel & Burrows 2012). The former, used by most investigators, predicts a higher luminosity for a given star at a given age with a given mass, than the latter. Neither set of models has been observationally confirmed for planetary mass objects at these young ages. Multi-M J companions, capable of driving observed two-arm spirals, may not be detectable by current surveys if the "cold start" models are appropriate -at 50 Myr, a 3M J hot start planet has roughly the same H-band luminosity as a 10M J cold start planet (Baraffe et al. 2003;Spiegel & Burrows 2012). Planets that are fainter than typically assumed would also help explain why direct imaging searches for the predicted perturbers in two-arm spiral disks have failed to detect massive giant planets at the flux levels predicted by the hot start models (SAO 206462 and MWC 758; Maire et al. 2017;Reggiani et al. 2018).
As an alternative way to resolve the discrepancy, one might imagine that multi-M J planets located at 10s to ∼100 AU from the star may be more common at a few Myr than at the 10-100 Myr ages typically probed by direct imaging studies (e.g., Bowler 2016), either because planets migrate inward from the original orbital distances via disk-planet interactions, and/or they are scattered out of the system via planet-planet interactions. Both options appear unlikely, however. While a planet could be driven inward by an outer disk located beyond the planet (e.g., Kley & Nelson 2012), in all systems with two-arm spirals the purported planet is located beyond the edge of the disk. In the planet-planet interaction scenario, multiple giant planets that are initially stabilized by the gas disk (e.g., Dunhill et al. 2013) later undergo planet-planet scattering that removes planets from the system once the gas disk dissipates (e.g., Dong & Dawson 2016). Systems with two-arm spirals, however, are unlikely to have additional giant planets with similar masses at distances close enough to interact with the purported arm-driving planet, as such planets will drive their own sets of two-arm spirals, which are not seen.
Thus, the most attractive explanation for two-arm spirals in disks is that they are driven by massive giant planetary companions ( 5M J ) that are fainter than predicted by hot start evolutionary models.
Planetary companions might also play a role in explaining multi-arm spirals. Numerical simulations have shown that a planet of mass ∼ M J usually drives one to three arms on one side of its orbit that are detectable in NIR imaging observations (Dong & Fung 2017). In -14disks with multi-arm spirals (Figure 1), the arms are at similar distances. To drive them with Jovian planets requires more than one planet, positioned over a small range of disk radii, potentially in conflict with theories of giant planet formation.

Summary and Future Prospects
Spiral arm structures are surprisingly common in NIR scattered light imaging of protoplanetary disks, especially among Herbig stars. Among the 11 Herbig stars that have been well-studied in scattered light, 6/11 have spiral arms of some kind, and 2/11 have near-symmetric two-arm spirals. Among the 24 single Herbig stars within 200 pc, 10 have been imaged in scattered light. The spiral arm occurrence rate in this volume-limited Herbig intermediate mass star (IMS) sample is ≥ 2/24 for two-arm spirals and ≥ 5/24 for spirals of all kinds; if some of the remaining 14 disks in the sample turn out to have spiral structures, the fractions will be higher. Although the sample studied to date is small and the occurrence rates of spiral arms poorly determined as a result, the data do have the potential to constrain quantities of fundamental importance for planet formation.
Two-arm spirals appear to be more common among Herbig IMS disks than are directly imaged giant planet companions to mature AB stars as estimated in the literature. This discrepancy is difficult to account for in a simple way given our current understanding of the mechanisms that can create spiral structure. Gravitational instability is unlikely to explain the two-arm spirals detected in class II disks mainly because of the short lifetime of this phase of disk evolution compared with the high fraction of two-arm spirals observed. The stellar accretion rates of Herbig systems with two-arm spirals are also not particularly high compared to other Herbigs. To explain the high occurrence rate of two-arm spirals, we propose that they may be produced by massive giant planets that are fainter than predicted by hot start models. Such planets must occur as companions to 8% of Herbig IMS.
Gravitational instability may drive multi-arm spirals if the disk mass is 10% of the stellar mass. Submillimeter continuum-based disk mass estimates for disks with multi-arm spirals are generally about one order of magnitude below this threshold. However, the high stellar accretion rates of Herbig stars, if sustained over a reasonable fraction of their lifetimes, suggest much larger disk masses. As a result, gravitational instability is a plausible explanation for multi-arm spirals. With the multi-arm phase of gravitational instability persisting longer than the two-arm phase, it is possibly consistent with the 3/11 occurrence rate among well-studied Herbig IMS and the ≥3/24 occurrence rate in the volume-limited Herbig IMS sample.
-15 -Several future studies can lend insights into the issues raised here. Our study is limited by the small sample size of disks that have been well-studied in scattered light. High signalto-noise scattered light imaging of the remainder of the volume-limited Herbig IMS sample would be extremely valuable in firming up the occurrence rates of spiral features. ALMA observations of the gas and dust emission structure of disks, which do not require illumination by the central star, are a valuable way to develop a more complete census of spiral arms. To distinguish between the various scenarios that could produce the observed spiral structures, it may be useful to compare not only the observed and expected number and symmetry of arms but also the observed and expected surface brightness contrast of spiral features above the rest of the disk.
The current sample is biased toward older sources that lack substantial infall and are therefore less likely to have gravitationally unstable disks. It would be very interesting to search for spiral structure among the younger precursor population (i.e., among GK-type intermediate mass T Tauri stars) either in scattered light or millimeter emission. The time evolution of the occurrence rate of two-arm spirals may shed light on the formation timescale of giant planets at large distances, and those of multi-arm spirals may reveal how disk masses evolve with time.
-36 -   , SAO 206462 (Garufi et al. 2013), LkHα 330 (Uyama et al. 2018, and DZ Cha . Multiple weak arms are observed on smaller scales in AB Aur (Hashimoto et al. 2011), HD 142527 (Avenhaus et al. 2017, and HD 100546 . A one-arm spiral is observed in V1247 Ori (Ohta et al. 2016). The scale bar in all panels is 0. 5. All images show polarized intensity, and are r 2 -scaled and on a linear stretch, except for HD 100546, which shows total intensity and is on a log stretch because the polarized intensity image does not reveal the spiral arms well ).   The disk mass is estimated from the submillimeter continuum, assuming a gas-to-dust mass ratio of 100. The source name corresponding to each number in the horizontal axis can be found in Table 2. All disk masses are < 0.2M , with the average M disk /M = 0.008. Herbigs with spiral structure are interspersed with other Herbigs. See §4 for details.    Sandell et al. (2011) performed simple model fitting to the SED of a sample of Herbig disks and found T dust ∼ 30−70 K. Andrews et al. (2013) pointed out the T dust dependence on the stellar luminosity L , and suggested scaling T dust as L 1/4 . Pascucci et al. (2016) suggested that the size of the dust disk should also be taken into consideration when estimating T dust .
Here, we introduce another correction factor on T dust based on the outer radius of the disk in submillimeter continuum observations, R mm (Column 12 in Table 2), as where the dust is emitting is also important -a dust ring further from the star is colder than a dust ring closer in, everything else being equal. We scale T dust as The factor of 0.5 in 0.5R mm is an approximate correction factor meant to identify the radius where the dust temperature is the global average. This is to account for the fact that dust at different radii has different temperature. If the surface density follows a Σ ∼ 1/R radial profile and the disk extends from R = 0 to R = R mm , 0.5R mm is the half-mass radius, i.e., the radius inside which half of the disk mass is enclosed. Equation B2 gives the specific temperature at R = 0.5R mm , not R mm , in a disk; this temperature is considered as a proxy for the average dust temperature in a disk whose outer edge is at R mm . For unresolved submillimeter sources, we adopted the average R mm in resolved sources: 191 AU for Herbig disks and 144 AU for T Tauri disks. This relation is calibrated against a model for HD 163296 , and reproduces the midplane temperature produced by several radiative transfer models of both T Tauri and Herbig disks in the literature to within 40% Dullemond & Dominik 2004;Dong 2015;Andrews et al. 2016).
Another simple scaling for estimating the disk mass comes from the stellar accretion rate (Hartmann et al. 1998). The stellar accretion rate declines roughly as t −3/2 ; thus the mass remaining in the disk is M disk = 2Ṁ t age . Because stellar accretion rates are time variable by a factor of ∼3, the mass estimate of an individual disk is uncertain by at least this factor. However, for an unbiased ensemble of disks, the mean mass of the disks is a reasonable estimate. For targets in our volume-limited Herbig IMS sample (Table 2), we conservatively estimate t age to be at least 2 Myr based on the HR diagram ( Figure 3).
-50 - Table 2 Measuring the accretion rate on to Herbig Ae/Be (HAeBe) stars is complicated by the fact that the stellar photosphere is bright at ultraviolet wavelengths. Thus the contrast between the accretion luminosity and stellar photosphere is low (e.g., Muzerolle et al. 2004). To address this challenge, several groups have calibrated the veiling of the Balmer jump against various emission line diagnostics for stars with the greatest ultraviolet contrast. In an X-Shooter survey of 91 Herbig Ae/Be stars, Fairlamb et al. (2015Fairlamb et al. ( , 2017 have calibrated the luminosity of 32 emission lines against the accretion luminosity as inferred from the veiling of the Balmer jump.

C. Accretion Rate of Herbig Stars in
We recalculated the accretion rates for all sources to obtain a self-consistent set of values. We adopted the accretion luminosity for our sources from the measurement of the Balmer discontinuity when this measurement was available as it is the most direct measure of stellar accretion. When the Balmer discontinuity was not available, we used the luminosity of Brγ or Hα (L Brγ , L Hα ) to calculate the accretion luminosity , Table 2): log L acc = (4.46 ± 0.23) + (1.30 ± 0.09) log L Brγ , log L acc = (2.09 ± 0.06) + (1.00 ± 0.05) log L Hα .
We prioritized Brγ over Hα because it is less sensitive to corrections for both reddening and photospheric absorption.
There are two caveats to keep in mind concerning the measurement of the stellar accretion rate of HAeBes. Firstly, it is not exactly clear how HAeBes accrete. While young A/B stars do not have the convective envelopes believed to play a key role in sustaining the kG magnetic fields observed among their later type counterparts, there is some evidence of magnetospheric accretion (Guimarães et al. 2006;Cauley & Johns-Krull 2014). If magnetospheric accretion controls the star/disk interface, it is likely through higher order field components. However, the conversion of the veiling of the Balmer discontinuity to an accretion luminosity relies on an assumption about the energy in the accretion column (F ) and the filling factor (f ; Muzerolle et al. 2004). Neither is well constrained directly for HAeBes. Adopting values that reproduce the observations for T Tauris stars may not be appropriate for these higher mass analogs. Secondly, early type stars with metal poor photospheres, in particular the λ Boo stars, show a UV excess (Murphy & Paunzen 2017). Folsom et al. (2012) and Kama et al. (2015) found a relatively high fraction of HAeBes in their samples (up to ∼50%) showing the λ Boo pattern to varying degrees. If this effect is common among the targets in our sample, it is possible thatṀ have been systematically overestimated. Future work to account for the uncertainty in F and f , as well as the effect of the photospheric abundance, will improve our understanding of the stellar accretion rate of these systems.

-51 -
The luminosities collected from the literature were rescaled using the distances available from Gaia Collaboration et al. (2016) when they were available. Conversion of equivalent width measurements of the HI lines to line luminosities requires correction for the veiled equivalent width of photospheric absorption, where W obs is the equivalent width of the emission line, W phot is the equivalent width of the photospheric line, W circ is the equivalent width of the circumstellar emission, and ∆m is the excess broadband emission above the stellar photosphere (Rodgers 2001). The equivalent widths of the photospheric features were taken from Fairlamb et al. (2017) when available. For stars not included in their study, stars with matching temperatures and gravities were used to infer the photospheric equivalent width. The veiling in the K-band was determined by the dereddened color excess, E(V − Ks). The intrinsic colors were taken from Pecaut & Mamajek (2013). For the sources where Hα was used, the extinction at 6563Å was inferred from E(B − V ), R V =5, and the reddening law presented in Mathis (1990).
The stellar accretion rate is proportional the accretion luminosity so that, where R is the stellar radius and M is the stellar mass. The stellar parameters were taken from the Siess et al. (2000) pre-main sequence evolutionary tracks and tabulated in Table 2.

D. Comment on Individual Systems
Stars for which the determination of the stellar parameters were not straight forward are discussed in detail. We also comment on the scattered light imaging observations of the AB Aur and HD 100546 disks.