Planet-driven spiral arms in protoplanetary disks: I. Formation mechanism

Protoplanetary disk simulations show that a single planet can excite more than one spiral arm, possibly explaining recent observations of multiple spiral arms in some systems. In this paper, we explain the mechanism by which a planet excites multiple spiral arms in a protoplanetary disk. Contrary to previous speculations, the formation of both primary and additional arms can be understood as a linear process when the planet mass is sufficiently small. A planet resonantly interacts with epicyclic oscillations in the disk, launching spiral wave modes around the Lindblad resonances. When a set of wave modes is in phase, they can constructively interfere with each other and create a spiral arm. More than one spiral arm can form because such constructive interference can occur for different sets of wave modes, with the exact number and launching position of spiral arms dependent on the planet mass as well as the disk temperature profile. Non-linear effects become increasingly important as the planet mass increases, resulting in spiral arms with stronger shocks and thus larger pitch angles. This is found in common for both primary and additional arms. When a planet has a sufficiently large mass ($\gtrsim$ 3 thermal masses for $(h/r)_p=0.1$), only two spiral arms form interior to its orbit. The wave modes that would form a tertiary arm for smaller mass planets merge with the primary arm. Improvements in our understanding of the formation of spiral arms can provide crucial insights into the origin of observed spiral arms in protoplanetary disks.


INTRODUCTION
A point mass perturber creates disturbances in the differentially rotating (e.g., Keplerian) background disk, which appear in the form of wakes spiraling away from the perturber. Such spiral structures are seen in observations and/or numerical simulations of a variety of astrophysical disks, including protoplanetary disks, circumplanetary disks, and disks around various binary systems (e.g., binary stars, binary black holes).
As a spiral arm steepens into a shock, its angular momentum is transferred to the background disk gas, opening a gap (Rafikov 2002b). When a single planet excites multiple spiral arms, it can create multiple gaps and pressure bumps in between . The resulting structure in the gas disk can produce corresponding, enhanced features in dust emission and scattering, such that a single planet can be responsible for more than one main gap. This mechanism possibly explains some of the multiple concentric gaps and rings seen in young protoplanetary disks (e.g., ALMA Partnership et al. 2015;Andrews et al. 2016;Tsukagoshi et al. 2016;van Boekel et al. 2017;Isella et al. 2016;Zhang et al. 2016;Ginski et al. 2016;Walsh et al. 2016). These structures can also affect subsequent planet formation in such ringed/gapped disks by collecting solid particles preferentially in local gas pressure maxima (i.e., rings) through aerodynamic drag.
Beyond the effects on dust and gas rings and gaps, spiral shocks can transport angular momentum and dissipate energy in various astrophysical disks including protoplanetary disks (Goodman & Rafikov 2001;Ju et al. 2016Ju et al. , 2017Rafikov 2016;Zhu et al. 2016;Ryan & MacFadyen 2017). Not surprisingly, this capa-bility is not limited to the primary spiral arm -the one directly attached to the planet -but the secondary spiral arm is also known to be able to contribute to angular momentum transport (Arzamasskiy & Rafikov 2017). If there exist additional spiral arms (e.g., tertiary, quaternary, ...) and they form through the same mechanism as the primary arm, angular momentum transport is also expected for the additional arms.
The characteristics of observed spiral arms can be used to constrain the masses of unseen planets. As we will show in a companion paper (Bae & Zhu 2017, hereafter Paper II) the number of spiral arms varies as a function of the planet mass. In addition, it is known that the arm-to-arm separation increases as a function of planet mass Fung & Dong 2015;Lee 2016, Paper II).
Despite the importance and growing number of numerical/observational studies showing multi-armed spirals in protoplanetary disks, the mechanism by which a planet excites multiple arms has not been fully understood. In Ogilvie & Lubow (2002, hereafter OL02), the formation of primary arm was explained as the result of constructive interference among a set of wave modes having different azimuthal wavenumbers. The gravitational potential of a point mass perturber Φ p at a position (r, φ) and time t can be decomposed into a Fourier series, a sum of individual azimuthal modes having azimuthal wavenumbers m = 0, 1, 2, ..., ∞: . Through the resonance between the rigid rotation of a perturbing potential and the epicyclic motions of disk material, the mth Fourier component of the potential launches m wave modes that are evenly spaced in azimuth (Goldreich & Tremaine 1978a,b, 1979, see also the review by Shu 2016). Throughout this paper, we use n = 0, 1, ..., m − 1 to represent each wave mode excited by the mth Fourier component. Among m wave modes launched by an arbitrary mth azimuthal component, OL02 considered the one that originates from the perturber position (hereafter n = 0 components; see Figure 1 for illustration). Using a linear wave theory, OL02 calculated the phases of n = 0 components having different azimuthal wavenumbers and confirmed with numerical simulations that these wave modes add coherently, creating a primary spiral arm.
In this paper, we show that additional spiral arms form in a similar way as the primary arm does: through constructive interference among appropriate sets of wave modes having different m. More specifically, in the inner disk, we show that the wave modes excited at the Lindblad resonance with an azimuthal shift of 2π/m from the n = 0 component (i.e., n = 1 components) create a secondary spiral arm. Similarly, a tertiary spiral arm forms via constructive interference among the wave modes excited at the Lindblad resonance with an azimuthal shift of 4π/m from the n = 0 component (i.e., n = 2 components), and so on and so forth when possible. In the outer disk, it is n = m − 1 components in each mth azimuthal component that form a secondary spiral arm. As we will show throughout this paper, this mechanism explains the characteristics of spiral arms known from previous numerical simulations very well. This paper is organized as follows. In Section 2, we first compute the phases of individual wave modes excited by a planet using a linear density wave theory and show that certain sets of wave modes can be in phase. We then carry out a suite of two-dimensional numerical simulations in Section 3 and verify that the sets of individual wave modes predicted by the linear theory add constructively on to each other, creating spiral arms. In Section 4, we highlight some important non-linear effects, based on numerical simulations with a range of planet mass covering three orders of magnitude. We summarize our findings and conclude in Section 5. In Paper II, we present a parameter study varying the disk temperature and planet mass and implications of the present work.

EXPECTATION FROM LINEAR WAVE THEORY
A rigidly rotating point mass perturber in a differentially rotating disk excites density waves at the vicinity of Lindblad resonances (Goldreich & Tremaine 1979). The perturbation driven by the mth azimuthal Fourier component of the perturber potential in a twodimensional (r, φ) plane can be written as whereX(r) is the amplitude of the perturbation, k(r) is the radial wave vector, Ω p is the orbital frequency of the perturber, and t represents time. The radial wavenumber k can be related to the azimuthal wavenumber m and the background disk properties through the WKBJ dispersion relation, which can be written in a two-dimensional disk as in the absence of self-gravity and dissipation processes. In the dispersion relation, Ω is the local orbital angular frequency, κ 2 ≡ (1/r 3 )d(r 2 Ω) 2 /dr is the square of the epicyclic frequency, and c s is the local sound speed. As a wave propagates, its phase varies over radius. The phase of a wave with an azimuthal wavenumber m at an arbitrary radius r, φ m (r), can be obtained by integrating dφ/dr = −k/m: where φ m (r 0 ) is the phase at r = r 0 . For a Keplerian disk (Ω ∝ r −3/2 , κ = Ω), one can re-write the dispersion relation in Equation (2) as follows: Inserting Equation (4) into Equation (3), with dφ/dr = −k/m, we obtain (5) Density waves driven by a point mass perturber launch around the Lindblad resonance r ± m = (1 ± 1/m) 2/3 r p , propagating inward in the inner disk and outward in the outer disk (Goldreich & Tremaine 1978b, 1979, where r p is the radius of the perturber's circular orbit. Far from the resonance, the phases of m wave modes are where n = 0, 1, ..., m − 1, as can be inferred from the asymptotic behavior of the Airy function (Ward 1986; OL02) 1 . Using Equation (5) and (6), we now obtain the phases of individual wave modes with any given combination of m and n: Throughout this paper, we call Equation (7) the phase equation.
We note that the phase equation consists of a constant component (the first two terms) that determines the launching position of wave modes in azimuth at the Lindblad resonance, and a radially varying component (the third term) that determines how tightly the wave modes are wrapped as they propagate away from the Lindblad resonance. The tightness of the wave modes depends on the azimuthal wavenumber m as well as the background disk temperature (c s ) and rotation profiles (Ω). The epicyclic term is often ignored in the literature, such that the dependence of the wave propagation on m is neglected. This may be a minor effect in many cases; however, the m-dependency in wave propagation is what enables the formation of multiple spiral arms by a single perturber, and therefore the epicyclic term has to be included.
To help visualize wave excitation and propagation, we present in Figure 1 a schematic diagram showing the phases of individual wave modes with m = 4 as an example.

Primary Spiral Arm Formation
1 One may use a phase offset term π/(3m) in Equation (6), which is the offset at exact Lindblad resonance locations (Ward 1986), instead of π/(4m). The difference between the two offset values are small (π/(12m)), especially when m ≫ 1 modes are considered, and we find that the formation mechanism of spiral arms is not affected by the choice of the offset value (i.e., π/(4m) vs. π/(3m)). In OL02, the formation of the primary spiral arm driven by a planet in a protoplanetary disk was explained as the result of constructive interference among wave modes having different azimuthal wavenumbers. The constructive interference considered was for n = 0 components of each mth azimuthal mode. Here, we follow OL02 and briefly summarize their findings since it will help understand the formation of additional spiral arms that will be explained in the following section. For the example presented in this section, we adopt a temperature profile that is decreasing as a function of radius following T ∝ r −1/2 (c s ∝ r −1/4 ). In addition, we limit our attention to a Keplerian rotation profile and adopt a disk aspect ratio at r = r p of (h/r) p = (c s /v φ ) p = 0.1 such that the sound speed is much smaller than the rotation speed.
For n = 0 components, φ m is independent of m in the large m limit as can be seen from the phase equation. This implies that waves with different azimuthal wavenumbers can have the same phase so constructive interference among the waves may be possible. In Figure  2, we present the phases of n = 0 components of wave modes having azimuthal wavenumbers m = 1−20, calculated with the phase equation. The phases are growing positively/negatively in the inner/outer disk, meaning that these wave modes are trailing waves. The fact that the phases become greater than 2π as they propagate to the inner disk, or smaller than −2π in the outer disk, indicates that these wave modes can wind up multiple times before they reach the disk inner/outer boundary.
In the right panels of Figure 2, we present the relative phases of n = 0 wave modes (φ m,0 ) to the phase of The upper panels present phases in the inner disk (r < rp) while the lower panels present phases in the outer disk (r > rp). Note that there is no m = 1 mode in the inner disk because its inner Lindblad resonance is located at r = 0. The 'I'-shaped marks in the right panels show the azimuthal width ∆φ ≈ 2π(h/r)p within which different modes have to be located to participate in the constructive interference. The dotted curves in the right panels show the phase of the primary arm in the numerical simulation with full perturber potential (see Section 3). m = ∞ wave mode (φ ∞ ) so that the phase difference among the wave modes can be more clearly seen. As can be seen from the figure, wave modes are nearly in phase when they launch; this is why the primary arm forms directly attached to the perturber. However, because wave modes with a small m are less tightly wound than the ones with a large m, as inferred from the phase equation, small m modes are left behind/ahead in the inner/outer disk as they propagate.
The perturbation driven by a point mass perturber in a disk is dominated by azimuthal wavenumber m ≈ (1/2)(h/r) −1 p (Goldreich & Tremaine 1980). In order for the wave modes with different m to be coherently added, they have to be within the wave crest generated by the dominating mode, ∆φ ≈ 2π(h/r) p . This azimuthal width is presented with 'I'-shaped symbols in Figure 2. As can be inferred from the figure, the constructive interference can fail for small m wave modes far from the planet.

Formation of Additional Spiral Arms
Extending the primary spiral arm formation scenario outlined in the previous section, we propose that the nth components of each azimuthal mode, where n is now non-zero, can become in phase as they propagate and form additional spiral arms through constructive interference. While n = 0 components launch almost in phase as seen in Figure 2, other n components launch with non-negligible phase differences. For example, if one would draw φ m,1 for different m in Figure 1, small m modes will launch with larger initial phases than large m modes. However, small m modes are less tightly wound than large m modes so it is possible that small m modes are caught up by large m modes as the wave modes propagate. As in Section 2.1 we compute the phases of different wave modes and examine whether or not constructive interference will be possible. We first focus on the inner disk in Section 2.2.1 and then move on to the outer disk in Section 2.2.2. Figure 3. Phases of n = 1 and n = 2 wave modes in the inner disk (r < rp), having different azimuthal wavenumbers from m = 1 (red) to 20 (purple): the left panels show the actual phase values while the right panels show the relative phase values to m = ∞ mode. Note that the shape of the phase curves for a given m are identical regardless of n, but only the launching point is shifted in azimuth for different ns. The 'I'-shaped marks in the right panels show the azimuthal width ∆φ ≈ 2π(h/r)p within which different modes have to be located for the constructive interference to occur. The dotted curves in the right panels show the phase of the secondary and tertiary arms in the numerical simulation with full perturber potential (see Section 3).

Inner Disk
In Figure 3, we plot the phases of n = 1 components for m = 2 − 20 azimuthal modes. As shown, small m modes launch at larger azimuthal angle, but large m modes catch up the small m modes in phase because small m modes are less tightly wound. In this specific example, m = 2 − 20 modes become in phase (∆φ 2π(h/r) p ) at r ∼ 0.3 r p . The same can happen for n = 2 components; however, small m modes in this case will launch at even larger initial azimuthal angles compared with n = 1 components, so the wave modes have to travel further in order to become in phase. For n = 2 components, m = 3 − 20 modes become in phase at r ∼ 0.1 r p . For the disk considered here, n > 2 components are unlikely to become in phase before they reach the disk inner boundary.

Outer Disk
We now turn our attention to the outer disk. We examine n = m − 1 and n = m − 2 components instead of n = 1 and 2 components, since the wave modes considered here are trailing waves and, again, small m modes are less tightly wound than large m modes. In Figure 4, we present the phases of n = m−1 and n = m−2 components. As seen in the figure, the phase differences among different m modes initially decrease, but remain nearly constant beyond r ∼ 3 r p . It hence appears that small m modes are not able to catch up to large m modes.
While we present the phases out to r = 5 r p only in Figure 4, constructive interference for n = m − 1 and n = m − 2 components beyond the radius is unlikely. This can be inferred from the phase equation. When r ≫ r p , the last term in the phase equation simplifies to and thus has no m dependence. This means that when different m modes launch at different azimuthal angles and they are not in phase before r ≫ r p , they will not be in phase in the outer disk. For the disk considered here, it is thus expected that only one arm forms in the outer disk through constructive interference among n = 0 components. On the other hand, more than one outer spiral arm can form when multiple sets of wave modes become in phase before r ≫ r p , which can occur in colder disks (see Paper II).

NUMERICAL SIMULATIONS: VERIFYING THE LINEAR THEORY PREDICTION
In Section 2, we showed that appropriate sets of wave modes having different azimuthal wavenumbers can be in phase from their launching points (n = 0 component) or as they propagate (non-zero nth components). In this section, we carry out numerical simulations to verify that constructive interference among the sets of wave modes predicted by the linear theory can indeed occur, generating spiral arms. We consider three models for this purpose.
Model 1: We carry out 20 calculations each of which includes a single mth Fourier-decomposed potential of a planet, where m = 1, 2, ..., 20. We then construct a single surface density output Σ by summing the perturbed density from each of the single mode calculation: where Σ init is the initial, unperturbed surface density and Σ m is the surface density obtained in a simulation with only the mth Fourier potential included.
Model 2: We carry out one calculation in which m = 1− 20 Fourier-decomposed azimuthal components of the planet potential are included.
Model 3: We carry out one calculation with the full planet potential.
By comparing Model 1 and 2 with Model 3, we will be able to test whether linear addition (i.e., superposition) of individual waves explains the main features of the model with the full potential. The comparison between Model 1 and Model 2 will allow us to examine if there exist any non-linear mode-mode interactions. If there is no interaction between different azimuthal modes at all, we expect that Model 1 and 2 will produce identical results. Finally, by comparing Model 2 with 3, we will be able to see the contribution from large m modes (m > 20) in generating spiral arms.  (7) is plotted with a black curve. Note the excellent agreement between the linear theory and numerical simulation.

Numerical Methods
We solve the hydrodynamic equations for mass and momentum conservation in the twodimensional polar coordinates (r, φ) using FARGO 3D (Benítez-Llambay & Masset 2016): In the above equations, Σ is the surface density, v is the velocity, P = Σc 2 s is the pressure where c s is the isothermal sound speed, Φ * = −GM * /r is the gravitational potential of the central star, and Φ p is the potential of the planet. The potential of the planet is where M p is the planet mass, r and r p are the radius vectors of the center of grid cells in question and of the planet, and s = 0.6 h p is the smoothing length. In this work, we ignore the indirect term which arises due to the offset between the central star and the origin of the coordinate system. The "full planet potential" in Equation (11) is used for Model 3. Assuming a circular planetary orbit, the potential in Equation (11) can be expanded into a Fourier series: Here, where δ ij is the Kronecker delta, β = r/r p , and b m 1/2 (β) is the Laplace coefficient defined as (Brouwer & Clemence 1961 The smoothing length s is included in the denominator of the right-hand side of Equation (14) in order for the summation of the Fourier-decomposed potential in Equation (12) to be consistent with the full potential in Equation (11). The Fourier-decomposed potential in Equation (12) is used for Model 1 and 2, with ms chosen following the model description.
We use the planet mass of M p = 0.01 M th , where M th ≡ c 3 s /ΩG = M * (h/r) 3 p is the so-called thermal mass (Lin & Papaloizou 1993;Goodman & Rafikov 2001), at which the Hill radius is comparable to the disk scale height. Assuming a solar-mass star, 0.01 M th is about 3 Earth masses. When M p ≪ M th , the excitation and the initial propagation of the density waves from the planet is known to be well approximated in the linear regime. Goodman & Rafikov (2001) predicts that spiral arms driven by a 0.01 M th planet steepen into shocks ∼ 6 scale heights away from the planet (see their equation 30). When we compare the phases of spiral arms driven by a 0.01 M th planet with the ones driven by a 0.001 M th planet in Paper II, we find that the difference in the phases of spiral arms are negligible not only within the ±6 scale height regions around the planet, but in the entire disk. This suggests that, although spiral arms driven by a 0.01 M th planet can steepen into shocks, non-linear effects are negligible. As shown in the following section, density waves excited by a 0.01 M th planet in numerical simulations indeed show an excellent agreement with the linear theory predictions.
Our initial disk has power-law surface density and   temperature distributions: Σ init (r) = Σ p (r/r p ) −1 and T (r) = T p (r/r p ) −1/2 , where Σ p and T p are the surface density and temperature at the location of the planet r = r p . We choose T p such that (h/r) p = 0.1, to be consistent with the disk model used in Section 2. The simulation domain extends from r in = 0.05 r p to r out = 5 r p in radius and from 0 to 2π in azimuth. We adopt 4096 logarithmically-spaced grid cells in the radial direction and 5580 uniformly-spaced grid cells in the azimuthal directions, with which ∆r : r∆φ ≃ 1 : 1. At the radial boundaries, we adopt a wave-damping zone (de Val-Borro et al. 2006) to suppress wave reflection. No kinematic viscosity is added in the simulations.

Simulation Results
We first present the perturbed density distributions δΣ m /Σ init , where δΣ m = Σ m − Σ init , from individual mode calculations (Model 1) in Figure 5. While Figure  5 includes results from m = 1 − 4 mode runs only, we note that the discussion below applies to all individual wave mode runs with m = 1 − 20.
Most importantly, the excitation and propagation of density waves in the numerical simulations show an excellent agreement with the linear theory. Each azimuthal component of the Fourier-decomposed potential excites m wave modes at the inner and outer Lindblad resonance; m = 1 mode does not excite any waves in the inner disk because the inner Lindblad resonance is located at r = 0. The perturbation from individual wave modes is < 1% over the entire simulation domain, supporting the fact that these waves are in a linear regime. The amplitude of the perturbation in all individual mode runs increases as the waves propagate, which is also in a good agreement with the expectation for linear waves (Rafikov 2002a).
In Figure 6, we display the perturbed surface density distributions δΣ/Σ init for Model 1, 2, and 3. The azimuthal distributions of δΣ/Σ init at r = 1.5, 0.6, 0.3, and 0.1 r p are presented in Figure 7 for more quantitative comparison among the models. All three models form three spiral arms in the inner disk and one spiral arm in the outer disk. The primary arm is directly attached to the planet, spiraling away from it. In the inner disk, the secondary and tertiary arms excite around r ∼ 0.3 r p and r ∼ 0.1 r p , respectively. Note that these radial positions are in a good agreement with the predictions made based on the phase argument in Section 2. The secondary and tertiary arms become narrower in azimuth and produce stronger perturbations as they propagate, indicating that the individual waves participating in the formation of these arms become closer in phase so the constructive interference become more effective. This is also consistent with the linear theory prediction (see e.g., Figure 3).
Comparing the models, we find that both Model 1 and 2 reproduce the full potential model (Model 3) fairly well. The major difference seen in Model 3 is that the primary arm is sharper and produces a larger perturbation close to the planet (e.g., r = 0.6 and 1.5 r p in Figure 7). This is because the Laplace coefficients (Equation (14)) that determine the strengths of the perturbation driven by individual azimuthal modes decline slowly with increasing m when β = r/r p is close to unity. The contribution from the azimuthal modes with m > 20 is therefore not negligible near r = r p . At the radius a spiral arm excites (e.g., 0.3 r p for the secondary and 0.1 r p for the tertiary) we see that all the three models agree with each other very well, suggesting that the excitation of additional spiral arms is a linear process. As spiral arms propagate, however, we see Model 2 and 3 deviate from Model 1. For example, the phase of primary arm in Model 2 and 3 is offset in phase from the primary arm phase in Model 1 at r = 0.1 and 0.3 r p . Also, at r = 0.1 r p the secondary arm breaks up into finer azimuthal scales in Model 2 and 3. This suggests that there could potentially be non-linear mode coupling (e.g., Lee 2016) even at this low level of perturbations.
To further ensure that it is the n = 0, n = 1, and n = 2 components from different azimuthal modes that generate the primary, secondary, and tertiary arms, we subtract each nth component from Model 1 one at a time when constructing the final surface density output. More specifically, we return the surface density associated with nth component in each individual azimuthal mode calculation to the unperturbed value as illustrated in Figure 8.
In Figure 9, we present the perturbed density distributions from Model 1 along with the ones from models without n = 0 components, n = 1 components (n = m − 1 components in the outer disk), and n = 2 components (n = m − 2 components in the outer disk). For more quantitative comparison among the models, we present the azimuthal distributions of the perturbed density at r = 0.1 r p in Figure 10. As seen in the figures the primary arm does not form when the n = 0 components are subtracted, the secondary arm does not form when the n = 1 components are subtracted, and the tertiary arm does not form when the n = 2 components are subtracted. Also, we note that removing the n = 0 components does not affect the secondary and tertiary arms. Likewise, removing the n = 1 or n = 2 components only affects the secondary or tertiary arms.
Previous studies pointed out that a negative density perturbation appears before the secondary spiral arm forms (e.g., Arzamasskiy & Rafikov 2017). We also find such a negative density perturbation in our simulations: as shown in Figure 6 and 7, the negative density perturbation just right side of the primary arm develops at r ∼ 0.6 r p and deepens inward before the secondary arm launches. In the constructive interference scenario we explain here it is obvious that, after n = 0 components (i.e,. wave crests) form the primary arm, the wave troughs between n = 0 and n = 1 wave crests have to be in phase before n = 1 wave crests become in phase to form the secondary arm. Similarly, the wave troughs between n = 1 and n = 2 wave crests become in phase before n = 2 wave crests form the tertiary arm, and this is what forms the negative density perturbation between the secondary and tertiary arms. In short, a negative density perturbation between spiral arms can be understood as constructive interference among wave troughs, as opposed to constructive interference among wave crests which forms a positive density perturbation (i.e., spiral arms).
Since the linear approach (i.e., superposition of individual wave modes) explains the formation and propagation of spiral arms well, we can make use of the linear wave theory to predict the phases of spiral arms. In the right panel of Figure 6, we present the phases of m = (1/2)(h/r) −1 p = 5 mode for n = 0, 1, and 2 components: φ 5,0 , φ 5,1 , and φ 5,2 . As shown, the phases of spiral arms predicted by the linear theory agree well with the phases of spiral arms in the numerical simulation.

NON-LINEAR EVOLUTION OF SPIRAL ARMS
As shown in the previous section, the linear wave theory explains the formation and propagation of spiral arms reasonably well for a sufficiently low-mass planet (i.e., 0.01 M th ). As the planet mass grows, however, non-linear effects are expected to play an increasingly important role. In order to investigate the non-linear effects, we run a set of simulations with various planet masses of M p = 0.1, 0.3, 1, 3 and 10 M th (0.1, 0.3, 1, 3, and 10 Jupiter mass assuming a solar-mass star), adopting the full planet potential as in Equation (11). All the other numerical setup except the planet mass remains the same as explained in Section 3.
In each of the simulations we introduce the planet at the beginning of the calculation with its full mass, instead of growing the planet mass over an extended period of time. We take this approach because planets with M p 1 M th open a gap around their orbit. The gap edges then become unstable to the growth of the Rossby wave instability (RWI; Lovelace et al. 1999;Li et al. 2000Li et al. , 2001, launching spiral waves that have different pattern speeds from the planet-driven spiral arms. For the planet masses considered here, we find that the RWI develops over about ten or more orbital times. By having the full planet mass from the beginning, spiral arms launched by the planet fully develop in the entire disk well before the RWI develops. This approach thus allows us to avoid the interference from RWI-driven spiral waves. With the background disk profile assumed here, the planet excites two or three spiral arms in the inner disk depending on its mass. In the outer disk, on the other hand, the planet excites only one spiral arm independently on the planet mass. To determine the phases of spiral arms from the simulations, we find the local maximum of the density perturbation in azimuth as we follow each spiral arm along radius. The phases of the primary, secondary, and tertiary arms for different planet masses are presented in Figure 11. In the figure, we also indicate the radial locations at which secondary and tertiary spiral arms start to shock disk gas. In order to diagnose the shock location, we compute the potential vorticity (PV) ζ ≡ (∇ × v)/Σ. The idea is that the PV experiences a jump at the shock front (Li et al. 2005;Dong et al. 2011;). In Figure 12, we plot the perturbed surface density distributions along azimuth at The two-dimensional distributions of the perturbed surface density δΣ/Σinit from Model 1. The primary, secondary, and tertiary arms are labeled with 'P', 'S', and 'T', respectively. In the other three panels, the distributions of δΣ/Σinit (left middle) without n = 0 components, (right middle) without n = 1 components at r < rp and n = m − 1 components at r > rp, and (right) without n = 2 components at r < rp and n = m − 2 components at r > rp are presented. Note that the primary arm does not form when the n = 0 components are subtracted, the secondary arm does not form when the n = 1 components are subtracted, and the tertiary arm does not form when the n = 2 components are subtracted. Figure 10. One-dimensional plots of δΣ/Σinit along azimuth at r = 0.1 rp. The black curves present δΣ/Σinit before the subtraction, while the red curves present δΣ/Σinit after the subtraction: (left) without the n = 0 components, (middle) without the n = 1 components, and (right) without the n = 2 components. Note that removing certain n components does not affect formation of other spiral arms. The primary, secondary, and tertiary arms are labeled with 'P', 'S', and 'T', respectively. Figure 11. The phases of (from left to right) the primary arm (φp) in the inner disk, secondary arm (φs), tertiary arm (φt), and primary arm in the outer disk, for planet masses of 0.01, 0.1, 0.3, 1, 3, and 10 M th . The dotted curves represent the phases for m = (1/2)(h/r) −1 p = 5 mode (φ5,n) calculated with the phase equation, while the dashed curves represent the phase for m = ∞ mode (φ∞). The arrows in the two middle panels show where each spiral arm starts to shock disk gas, diagnosed based on the potential vorticity jump. In the middle right panel, the yellow dashed curve presents the primary arm phase for the Mp = 3 M th case (see text). The distortion seen in the outer primary spiral arm for 10 M th planet is due to too strong shocks at the arm front.
some selected radii to show the level of perturbation driven by spiral arms and the morphology of spiral arm front. We highlight some important aspects of spiral arm formation and propagation below.
First, the phases of secondary and tertiary arms during their excitation and initial propagation agree reasonably well with the linear wave theory prediction. Looking at the phases of the secondary arms in Figure 11 first, one can see that they are very closely located to each other in azimuth at r = 0.4 − 0.6 r p for such a broad range of planet mass, and moreover follow the linear prediction (φ 5,1 ) very well. This is also clear in Figure 12: at r = 0.4 r p , the secondary spiral arms are located close to the linear theory prediction. Note also that secondary arms from planets with masses of ≥ 0.3 M th have already evolved into shocks at this radius, as the steep density gradient as well as shock locations presented in Figure 11 suggest. The fact that spiral arms follow the linear theory well after they start to shock disk gas supports that shocks need to propagate some distance before they deviate from the linear theory (Goodman & Rafikov 2001). Similarly, the tertiary arms are closely located in azimuth at r = 0.15 − 0.2 r p and the linear theory predicts the phases of the tertiary arms very well at the radii.
Second, spiral arms are more opened with a larger planet mass. Spiral arms deviate from the linear theory prediction after their initial propagation, which ends sooner (i.e., at larger radii in the inner disk and at smaller radii in the outer disk) for larger planet masses (Goodman & Rafikov 2001; see also Figure 11). When a spiral arm non-linearly steepens into a shock it travels at a faster speed, resulting in a less tightly-wound shape than the linear theory prediction. The speed of the shock expansion is proportional to the amplitude of the shock, so spiral arms excited by a more massive planet propagate at faster speeds and thus appear to be more opened . In Figure 13, we present the measured pitch angle of spiral arms driven by 0.1, 1, 10 M th planets from numerical simulations, along with linear theory predictions. As expected, the pitch angles measured in simulations with 1 and 10 M th planets are larger than the linear theory prediction, while the linear theory and simulation agree well with each other for a 0.1 M th planet. Compared at a given radius, the more massive the planet is, the more opened a spiral arm is in general. This trend is commonly seen, not only for the primary arm but also for additional arms. We confirm this trend in disks with other (h/r) p values from the parameter study carried out in Paper II (see their Figure  5).
Third, only two spiral arms form in the inner disk for sufficiently large planet masses ( 3 M th ). We conjecture that this is because the primary arm non-linearly propagates and merges with the wave modes that would form a tertiary arm for smaller mass planets. In the third panel of Figure 11, we over-plot the phase of the primary arm driven by a 3 M th planet. As shown, the For a visualization purpose, each curve is scaled by a factor presented on the right side of each panel such that the peak of the primary arm at the radius has δΣ/Σinit = 1. The primary, secondary, tertiary arms are indicated with 'P', 'S', and 'T', respectively. The dotted vertical lines present φ5,0, φ5,1, and φ5,2 at each radius, while the dashed vertical line presents the phase of m = ∞ mode φ∞ at the radius.
primary arm becomes more opened as it non-linearly propagates to smaller radii and eventually overlaps with the linear phase of the tertiary arm at r ∼ 0.2 r p . At r = 0.3 r p in Figure 12, the primary arm and the tertiary arm are well separated in azimuth for a 1 M th planet. For a 3 M th planet, on the other hand, the primary arm has a broad wing-like density enhancement on the left side of the 'N'-shaped shock, which is likely formed by the waves that would constructively interfere with each other to form a tertiary arm at the azimuth. The fact that the primary arm produces a comparable magnitude of perturbation to the secondary arm with M p = 3 and 10 M th at r = 0.2 and 0.3 r p , whereas the secondary arm is generally stronger for lower-mass planets at these radii because the primary arm weakens due to less efficient constructive interference, also supports the idea of tertiary arm-forming waves merging with the primary arm.
Last, interference between spiral arms may occur. In Figure 14 we present the two-dimensional density distribution from the M p = 1 M th model. Also presented with a black curve in the figure is the predicted primary spiral arm front position based on a non-linear shock expansion theory (Goodman & Rafikov 2001;Rafikov 2002a), following the procedure detailed in Section 4 of Zhu et al. (2015). Note that the shock expansion theory predicts the primary arm phase reasonably well at r 0.25 r p , but fails inward of the radius where the primary arm gradually approaches to the secondary arm in azimuth. The primary arm then returns back to the predicted position at r 0.08 r p . Interestingly, we find that this is not a transient but a long-lasting and stationary feature. We propose that one possibility for this is the interference between the primary and secondary spiral arms. Figure 14. The two-dimensional distribution of the perturbed surface density δΣ/Σinit from Mp = 1M th model. The black curve presents the predicted shock front obtained with a non-linear shock expansion theory (see text). The black and red arrows on the right indicate the primary and secondary arms, respectively. Note that the primary arm deviates from the non-linear shock expansion theory prediction at ∼ 0.25 rp approaching to the secondary, and then returns back to the predicted phase at ∼ 0.08 rp.
The perturbation driven by the primary arm is expected to gradually decrease since it dissipates while propagating. However, the primary arm gains in strength inward of ∼ 0.2 r p as it approaches to the secondary arm while the secondary arm loses its strength over the same radii, possibly because some low azimuthal modes that become out of phase from the secondary arm are added to the primary arm and/or because of non-linear mode-mode interaction. This supports the hypothesis that the deviation of the primary arm from the shock expansion theory at 0.08 r p r 0.25 r p is due to the interference from the secondary arm. In Paper II, we find that interference between spiral arms in colder disks can even result in merging of the arms, presumably made possible because spiral arms in colder disks launch with smaller azimuthal separations. While we see potential evidence of interference between spiral arms, it is unclear at the moment under which conditions such interference occurs. Interference between spiral arms is an interesting phenomenon to study, but it is beyond the scope of the paper and thus a more thorough investigation is deferred to a future paper.
One thing that does not change regardless of planet mass (and also disk temperature; see Paper II) is that small m modes are less tightly wound than large m modes. This property of waves suggests that an additional spiral arm always forms ahead of the previous arm in azimuth in the inner disk and behind the previous arm in azimuth in the outer disk. Also, an additional spiral arm always forms farther away from the planet in radius, because the wave modes have to travel a longer distance to be in phase.

SUMMARY AND CONCLUSION
We have shown how a planet excites multiple spiral arms in the underlying protoplanetary disk. Using a linear wave theory we first calculated the phases of individual wave modes excited by a planet and showed that appropriate sets of wave modes having different azimuthal wavenumbers can be in phase, from their launching points (in case of the primary arm) or as they propagate (in case of additional arms). By carrying out a suite of two-dimensional hydrodynamic simulations, we then verified that the sets of wave modes predicted by the linear theory add constructively on to each other and form spiral arms.
As the planet mass grows, non-linear effects play an increasingly important role. We investigated the nonlinear effects by carrying out numerical simulations with various planet masses from 1 % of a thermal mass to 10 thermal masses.
Our main findings are: 1. The formation of spiral arms -both primary and additional arms -is a linear process: constructive interference among appropriate sets of wave modes having different azimuthal wavenumbers.

2.
A planet excites m evenly spaced wave modes at the mth Lindblad resonance, where m = 1, 2, 3, ..., ∞ is the azimuthal wavenumber. Among the wave modes the n = 0 components, which are the ones originating from the planet location, add constructively on to each other and form the primary arm (Section 2 and 3), confirming the mechanism presented in OL02.
3. Additional spiral arms form in a similar manner to the primary arm, but through constructive interference among non-zero nth components. Nonzero nth components excite out of phase at the Lindblad resonance, in contrast to n = 0 components, but constructive interference among the wave modes is possible because wave propagation is dependent upon their azimuthal wavenumber (Equation 7).
4. Phases of spiral arms follow the dominating azimuthal mode with m ≈ (1/2)(h/r) −1 p reasonably well, until they steepen into shocks and depart from the linear regime. We provide a generalized analytic formula in Equation (7), which can be used for the primary arm but also additional arms.
5. In the outer disk, the propagation of wave modes becomes independent on azimuthal wavenumber m when r ≫ r p . Additional spiral arms can thus form in the outer disk only if wave modes become in phase before r ∼ a few × r p (Section 2.2.2).
6. Spiral arms excited by a more massive planet propagate at faster speeds and thus appear to be more opened, consistent with previous studies (e.g., Zhu et al. 2015). This trend is seen in common, not only for the primary arm but also for additional arms (Figure 13).
7. Only two spiral arms form for sufficiently large planet masses (M p 3 M th ). The wave modes that would form a tertiary arm for smaller mass planets merge with the primary arm ( Figure 11 and 12).
To conclude, the multiple spiral arm formation mechanism presented in this paper can explain many characteristics of planet-driven spiral arms known from previous studies. An obvious extension is to examine this scenario in three dimensions, particularly when the background disk is vertically stratified. Because of vertical gravity and/or buoyancy, wave modes in such disks will behave differently as they depart from the disk midplane. Spiral arms thus may not have a coherent vertical structure, possibly explaining the curvature of spiral arms seen in three-dimensional numerical simulations (e.g., Zhu et al. 2015;Bae et al. 2016b). While we focused on the case in which the primary object (i.e., star) is much more massive than the perturbing companion (i.e., planet), the spiral arm formation mechanism presented here can also be applied to the systems where the companion body has a comparable mass to the primary (e.g., disks around dwarf novae or binary black holes) or a much larger mass than the primary (e.g., circumplanetary disks). In the cases when the companion body has a mass comparable to or greater than the primary mass, it is very likely that the companion launches largely open two-armed spirals in the disk around the primary. We present applications of the present work in Paper II and discuss whether various characteristics of observed spiral arms can be used to constrain the masses of yet unseen planets and their positions within their disks.