The Kepler-223 resonance holds information on turbulence during the gas disk phase

Many fundamental physical processes regarding planetary formation in protoplanetary disks are still imperfectly understood, with an elusive phenomenon being turbulence in such disks. Observations are available of planetary systems and of some protoplanetary disks, which can serve as a starting point for these investigations. Detected systems reveal different architectures of planets. One particularly interesting case to consider is the Kepler-223 system, which contains a rare configuration of four planets in a resonance chain, implying a certain migration history. We aim to use the orbital configuration of Kepler-223's planets to constrain the parameters of the protoplanetary disk that allow for the formation of a chain of mean-motion resonances that resembles the one of Kepler-223. The parameters we aim to investigate are primarily the disk viscosity and surface density. We use the swift_symba N-body integrator with additional dissipative forces to mimic planet-disk interactions. We constrain the surface densities and viscosities that allow the formation of a resonant chain like that of Kepler-223. We find that surface densities of up to a few minimum mass solar nebula (MMSN) surface densities and disk viscosity parameters $\alpha$ of a few $10^{-3}$ up to $10^{-2}$ are most successful at reproducing the architecture of this particular planetary system. We describe how these two quantities are linked to each other considering the success of reproducing the chain and find that higher disk surface densities in turn require lower viscosities to build the chain. Our results show that well characterized observed planetary systems hold information about their formation conditions in the protoplanetary disks and that it is possible to extract this information, namely the initial disk surface density and viscosity, helping to constrain planet formation.


Introduction
The study of exoplanetary systems is a key element for understanding the formation and evolution of planets, as well as the assembly of planetary systems and their dynamical evolution. Currently, the exoplanet sample is largely dominated by data from the Kepler mission, which provided well constrained planetary radii and orbital periods using the transit detection method. Out of all detected planets, small and close-in planets, the so-called Super-Earths and Mini-Neptunes, are thought to be orbiting up to 50% of sun-like stars (Howard et al. 2012;Mayor et al. 2011;Mulders et al. 2018), and they often occur in multiplanetary configurations, where the planets within each system mostly have similar masses and radii (e.g., Millholland et al. 2017;Weiss et al. 2018).
A notable example is the planetary system hosted by Kepler-223, a G5V star with a mass of 1.1M . Four Sub-Neptunian planets in mean motion resonance (MMR) are orbiting it (Mills et al. 2016). In addition to orbital periods and radii, Mills et al. (2016) used transit timing variations (TTVs) to obtain masses, eccentricities and inclinations of the planets. The masses are given by 7.4 +1.3 −1.1 M Earth , 5.1 +1.7 −1.1 M Earth , 8.0 +1.5 −1.3 M Earth and 4.8 +1.4 −1.2 M Earth , for planets b-e, respectively. The period ratios of the planets are 8:6:4:3 starting from the inner-most planet Kepler-223 b. This suggests the existence of individual MMRs of the pairs of planets in the MMR chain, that is a 4:3 resonance for the first pair, a 3:2 resonance for the second pair and a 4:3 resonance for the third pair. Mills et al. (2016) confirmed the resonant nature of the system, showing librating Laplace angles. More detailed information about the planets orbiting Kepler-223 is presented in Table 1.
While capturing into a MMR is a natural outcome of type-I migration, they are found very infrequently in the exoplanetary sample. It was suggested that capture into a resonance chain does not produce dynamically stable systems often, if the disk in inviscid (McNally et al. 2019). Another often studied reason for the rarity of MMRs is their instability during (e.g., Liu et al. 2017) or after the dispersal of the protoplanetary disk (e.g., Izidoro et al. 2017Izidoro et al. , 2019Pichierri & Morbidelli 2020;Lambrechts et al. 2019). Thus, the existence of the particular configuration of the Kepler-223 system, as well as other resonant chains like Trappist-1 (Gillon et al. 2016(Gillon et al. , 2017Luger et al. 2017) or GJ-876 (Marcy et al. 2001;Rivera et al. 2005), poses a definite challenge to formation scenarios of planetary systems. For Kepler-223, the resonant state appears to be long-term stable (it Article number, page 1 of 12 arXiv:2110.13835v1 [astro-ph.EP] 26 Oct 2021 A&A proofs: manuscript no. main Table 1: Parameters of the planets orbiting Kepler-223 (Mills et al. 2016 is, for example, protected from instabilities due to the emergence from secondary resonances (Pichierri & Morbidelli 2020)), raising the question of what condition brought about this particular configuration. Indeed, these exotic orbital configurations offer an unique opportunity, as the current state can be used to constrain parameters of models that describe the protoplanetary disk in which they formed. Important disk parameters that could be constrained this way are, for example, the surface density or a description of turbulent viscosity in the disk. The surface density is a measure for the mass, which in turn is an important parameter for the migration speed (Kley & Nelson 2012;Baruteau et al. 2014), and also relevant for example for gravitational instabilities. Turbulence in disks can arise from different hydrodynamical instabilities, for instance the global baroclinic instability (Klahr & Bodenheimer 2003) or the vertical shear instability (Richard et al. 2016;Flock et al. 2017;Manger et al. 2020), or from the magneto-rotational instability (MRI) (Balbus & Hawley 1991). Turbulent viscosity is a key element of many planet formation scenarios and disk evolution models, which is solidified by the Atacama Large Millimeter/submillimeter Array (ALMA) data of disks that might show still ongoing planetary formation Pinte et al. 2018;Keppler et al. 2018), and is typically described in terms of the α parameter (Shakura & Sunyaev 1972). Despite its significance and recent improvements, the α parameter is not very well constrained observationally, with values typically ranging between 10 −4 and about 4 × 10 −2 , as studies of disks observable by ALMA have shown (Rafikov 2017;Dullemond et al. 2018;Flaherty et al. 2017Flaherty et al. , 2018. These values are also consistent with hydrodynamical and MRI simulations (e.g., Turner et al. 2014).
In this work, the fact that turbulence causes stochastic fluctuations of the disks density structure, leading to stochastic forces on the planets, is used to study its effects on the formation of Kepler-223's resonance chain. The main idea is that turbulence can be a reason for destroying or preventing the formation of a MMR chain even before the disk has dissipated (e.g., Batygin & Adams 2017). The aim of this investigation is to constrain the α parameter in combination with other disk parameters, just from observational information contained in the current state of the planetary system.
The methods used to investigate Kepler-223 are laid out in Sect. 2. Kepler-223 serves as a test case for the methodology that we use to find constraints on the disk turbulence and surface density. However, the methods we apply are general and can be applied to other systems with resonant chains. The results of our investigation are presented in Sect. 3. We discuss implications and shortcomings of our investigation in Sect. 4.

Methods
We reproduced the Kepler-223 resonance chain using N-body simulations using the swift_simba code, with additional dissipative forces implemented to mimic the interactions of the plan-ets with the protoplanetary disk, including the effects of stochastic surface density changes caused by turbulence.

Type-I migration
As a model for the effects on the semi-major axis and the eccentricity resulting from planet-disk interactions, the prescription developed by Cresswell & Nelson (2008) was implemented. The effects of eccentricity damping and type-I migration are viewed as two separate contributions. The former is described bẏ with the timescale τ e given by Here, τ wave is the typical type-I migration timescale, where q = m M is the fraction of the planetary mass to the stellar mass, h = H/r is the disk's aspect ratio and the subscript p denotes that the quantity is evaluated at the planets location. Planetary type-I migration is modeled given the negative torque, with a migration timescale of again in the limit of small eccentricities, with the power law of the surface density α Σ (see below). The resulting change of the semi-major axis can then be calculated aṡ with τ a = τ mig /2 and p ≈ 2 for small e. Note that the migration speed scales linearly with q and the surface density Σ. The ratio of the migration to eccentricity damping time-scale, the K-factor (Ramos et al. 2017) is given by Equation (7) shows that, for typical values for h and α Σ , eccentricity damping happens on a much shorter timescale than migration. Therefore, in this regime the planets are expected to have vanishing eccentricities before capturing into resonance.
To stop migration at the inner edge of the protoplanetary disk, a planetary trap is implemented. To model the sign flip of the torque and its radial dependence, which ultimately stops planetary migration once the planet reaches the inner edge (Masset et al. 2006;Flock et al. 2019), a multiplier F is applied to the acceleration of the planets (Pichierri et al. 2018). It is given by Here, r e is the location of the edge (in AU) and h e is the aspect ratio of the disk at r e . This implementation does not affect eccentricity damping.

Disk-driven turbulence
Our prescription for disk turbulence follows an implementation of Ogihara et al. (2007) based on results from hydrodynamical simulations investigating MRI turbulence (Laughlin et al. 2004). Turbulent modes are generated with a turbulent strength parameter γ. We summarize this method in the following. The stochastic force exerted by the turbulent disk is described by where the parameter Γ is given by The potential is described by a superposition of individual modes, In equation (12), ξ describes a Gaussian distribution with σ = 1 and (r, θ) the radial and azimuthal position of the planet, while r c and φ c describe the center of the density fluctuations caused by the turbulent modes. The azimuthal extent of the modes is given by 2πr c m , while the radial extent is given by σ = πr c 4m . The time dependence of the modes is given by ∆t. A mode is created at time t 0 and dies when ∆t >t t − t 0 . For ∆t, the sound crossing time in angular direction (∆t = 2πr c mc s ) was chosen. To create a new mode, r c is randomly generated inside the considered edges of our system r in and r out and φ c ∈ [0, 2π) is also chosen randomly. The wavenumber m is chosen from a logarithmic random distribution in the range of 2 ≤ m ≤ 64. The resulting torque Γ t is then given by where Λ s,m = ξe Here, m pl denotes the planetary mass in order to avoid confusion with the wavenumber m.
The stochastic force strength parameter γ can be related to the α parameter typically used to describe turbulence in current disk models. To this end, we use on the one hand a relationship for the diffusion coefficient in the eccentricity D e to the model's strength parameter (Ida et al. 2008;Okuzumi & Ormel 2013), where the numerical pre-factor used was modified to fit the setup used here (see appendix A), on the other hand we use a similar relationship involving the α parameter (Okuzumi & Ormel 2013), This yields a conversion between the two parameters that reads This relationship should be considered valid as an order-ofmagnitude correspondence between γ (the parameter used in the disk model for generating turbulent waves) and α (the commonly used turbulence parameter).

Treatment of individual planet pairs
As a starting point, we treated all three planet pairs of the Kepler-223 system separately. We first simulate each resonant pair separately for varying surface density values Σ = Σ 0 (r/r 0 ) −α Σ in order to get a first overview over how the surface density distribution affects each individual pair. For all conducted simulations here and in later sections, the aspect ratio was fixed at h = 0.05. The power law of the surface density α Σ was kept constant at a value of α Σ = 1.5. The reference distance r 0 was chosen to be r 0 = 1 AU. For the purposes of this step, we set the trap at the inner disk edge at the position of the inner planet in each pair under investigation. This allows the inner planet to be fixed allowing a clear investigation of the convergent migration behavior, and in turn means that Σ 0 sufficiently describes the surface density. In fact, the planets also capture in order, from the inner-most to the outer-most one, when considering the case of all four planets, so this is a good simplification to make to gain first insights into the migration behavior. Without considering turbulence and for fixed planet parameters, the migration timescale is the quantity that mainly defines the migration behavior and speed, and thus the final resonance. It is a function of the surface density, which makes it an important quantity to investigate for migration. It is important to keep in mind that we do not actually know at which distance from the star the four Kepler-223 planets formed, nor when. Therefore, for the purposes of this initial step, the planets of each pair were put wide of the 2:1 MMR initially, only in order to find a rough first estimate for the values for Σ 0 that we expect to work for the formation of the specific resonance chain. As a general rule, for a higher surface density, resonances with lower k are skipped, subsequently leading to a capture into a resonance with a higher k (Batygin 2015). Therefore, we consider a range of values for In red, the error is indicated, which is related to the mass error of the planets (Mills et al. 2016). The black dotted lines are shown to make clear that within the errors, no value for Σ 0 can be found that allow each separate pair to capture into the correct resonance, if the pairs are initially wide of the 2:1 MMR. The blue dotted vertical lines indicate that this does not hold if the planets are instead placed just wide of the correct MMR, which eliminates the lower limit for Σ 0 .
the pairs capture into. Considering the different location and different masses of the pairs, we expect a value for Σ 0 to lead to different final k for each individual pair. We then compare these final k to the resonance indices measured by Mills et al. (2016) to find a range of surface densities that allow capture into the correct resonance. The result of this consideration is shown in Fig. 1. Under the strong assumption that the planets formed outside of the 2:1 MMR, one cannot find a surface density power-law prescription that allows all pairs to individually capture into the correct resonance. While in principle other slope-values α Σ than the used value of 3/2 are possible, this statement is independent of this parameter. However, when placing the individual pairs in front of the desired resonances, there is no lower limit for the surface density that allows the planets to be captured into that resonance. In that case, one can indeed find surface density values that allow the formation of the observed resonance for all pairs individually. All in all, this investigation is of course insufficient for the problem at hand, and should rather be taken as an estimate for the relevant range of Σ 0 for a more complete consideration. Naturally when considering a resonant chain of four planets, fourplanet dynamics need to taken into account, which also change the outcome with respect to the final configuration per surface density.

Treatment of all four planets together
The next natural step to investigate the capture into a resonant chain of the Kepler-223 system is thus to consider all four planets at once. We consider different surface densities, expecting it to influence the migration timescale, and different initial planet locations. For the planets' locations, we take a pragmatic ap-proach. As we already mentioned their initial formation location can not be inferred from observations and could have been in principle in a configuration where individual planet pairs lie far outside the currently observed resonances. For the purpose of this investigation, however, the planets are placed just outside the resonances in which they are observed today. The reason for this is that, since the innermost planet is locked at the inner edge of the disk (Masset et al. 2006;Ataiee & Kley 2021), migration is typically convergent anyway (see below), and thus if a planet pair is observed today in a certain mean motion resonance, it is natural to assume, within the framework of our investigation, that they once lied just wider of that specific period ratio. Note that this initial placement does not pose a contradiction to the results of Sect. 3.1.1. In that section, the premise was to find a range of potential surface density values that locally allow the formation of a resonance with index k that matches the observations by Mills et al. (2016). This was achieved by considering the analytically well-understood case of a two-planet system without turbulence. In the investigation in this section, the dynamics of the four-planet system while including effects of turbulence are considered. In this scenario, planets are able to reach the position where we put them at the start of the simulation, even if they do form wider of the resonances they are presently observed in and have to migrate inward. The distance to an unperturbed semimajor axis ratio is parameterized by ε, where i = 1, 2, 3 denotes the planet pair of the chain, and k i is the resonance index, with k 1 = 4, k 2 = 3 and k 3 = 4. The first planet's initial semi-major axis is set to a 1 = 0.076 AU, which is close to the current location of Kepler-223b given by Kepler's third law, a b = 0.0772 AU (Mills et al. 2016).
We run N-body simulations for values of Σ 0 ranging between 0.1 and 10 Σ MMSN , for all combinations of values for ε i out of {0.02,0.035,0.05} representing small, medium and large initial separations between the planets. Even for the largest initial separation, the planet pair was always placed in such a way that no first-order MMR had to be crossed in order to reach the intended one. After the integration, a system was considered to exhibit the correct chain of resonances if their period ratios were wide of the observed resonance, as well as the corresponding resonant angles showing libration. A system's configuration of periods was considered to be compatible with observations if the ratio of the periods of each individual planet pair exceeded the nominal value, but with a relative error below 0.1%. The result of this investigation can be seen in Fig. 2. One can see the surface density value at which the migration timescale is too short to be able to form the Kepler-223 resonance chain, because the needed resonances are skipped as expected. This happens beyond a surface density corresponding to ∼6 Σ MMSN . This demonstrates at face value that, at least without any treatment of turbulence, the minimum-mass extrasolar nebula (MMEN, Chiang & Laughlin 2013), cannot accurately describe protoplanetary disks of extrasolar systems like Kepler-223. It predicts a surface density of ∼7Σ MMSN at 0.1AU, which would not allow the formation of the Kepler-223 system.
From the figure, it appears like the capturing process of this four-planet system into resonance is chaotic and therefore depends very sensitively on both the initial separation and the surface density, where even a small change in the initial distance to the resonance may make the difference between a successful formation of or the failure to create the correct chain. Configurations where the correct chain was built are indicated by green triangles, while those who failed to do so are shown as red crosses. The vales for ε are ordered such that ε 1 is the top-most value. Simulations whose symbol is encircled share their initial parameters with those considered for a consideration including disk turbulence (Sect. 3.2).

Results with turbulence
Some of the initial conditions that lead to the successful creation of the Kepler-223 chain are used to investigate the effects of turbulence on the formation of the chain. Three different initial separations were picked, where surface densities ranged between 0.72 Σ MMSN and 3.93 Σ MMSN . These surface densities were chosen because they allow for most of the successful simulations. Choosing higher surface densities give too high of a migration speed, which is unfavorable for the formation of the resonance chain, while lower surface densities do not offer great constraints for the simulations that include turbulence. The initial conditions that were used are denoted with a blue circle in Fig. 2.
After finding initial conditions that are favorable for the formation of the resonance chain, we investigated the effects of turbulence on its formation to find constraints for the turbulence strength present in the protoplanetary disk during the formation of the chain. We based our analysis on 100 simulations per initial setup. The integration time is not based on physical properties of the system, since the remaining disk lifetime at the point of formation of the chain is not known, and neither is the influence of a possible time variability of α. Rather, the integration time was chosen as a compromise between simulation run time and stability of the results; we found that integrating for 10 5 years achieves this compromise. We also investigated the effects of a longer integration for a particular case, and found that, while some initially stable configuration became unstable, the difference was below 10%. Therefore, it does not significantly change the results presented in Fig. 3.
Four different strengths of turbulence were considered, ranging between α = 10 −4 and α = 5 × 10 −2 . The resulting trend is depicted in Fig. 3. Several trends are worth pointing out. Firstly, a larger turbulence strength generally leads to less success in forming the Kepler-223 resonance chain. The point at which the formation of the chain becomes significantly improbable lies at lower α for higher surface densities, which fulfills the expectation that the strength of the random kicks be dependent on the product Σ √ α (cfr. Eqs. (13), (17)). This can be further seen by considering individual columns of the figure, where a larger amount of simulations lead to the formation of the resonance chain when the surface density is lower at constant α. However, this trend does not hold for the lower turbulence strengths. This can be seen as an indication that simple analytical trends derived for the two-planet case (Batygin & Adams 2017) are not applicable when the strength of the turbulent kicks is not suffi-cient to disrupt chains that have already established their resonant configuration (see below), while they can serve as an orderof-magnitude approximation in a case where the turbulence is strong enough to disrupt formed chains. In the former case, disruption has to take place during the capturing process, where simple assumptions do not hold for four planets.
Systems can fail to recreate the Kepler-223 resonance chain for a variety of reasons. Frequent examples include the second planet pair skipping the desired 3:2 resonance and capturing into a 4:3 resonance instead, which can lead to other pairs failing to capture into or stay in the correct resonance. The second planet pair is especially prone to temporary capture, because its desired resonance has the lowest index k and its outer planet is more massive than its inner one (Deck & Batygin 2015). Even in the event of successful capture of the second pair, planets can skip their desired resonance during the initial capturing process of the chain or break out of it later on. Another possible path to failure is created by the brief divergent migration of the third pair, which occurs because its inner planet is more massive than the outer one. If the more massive planet is not captured fast enough into the 3:2 resonance of pair 2, pair 3 can cross the 3:2 resonance from the inside, which is known to cause a rapid excitation of eccentricity. The rapid excitation usually breaks the resonances of other pairs. The frequency at which the third pair crosses the 3:2 resonance is related to the initial separation of the planets. Once a desired resonance is skipped, the planets may either capture into a resonance with higher index k, or collide. Rarely, a pair captures into a second-order resonance. Figure 3 also contains the results of a simulation corresponding to initial conditions that, without considering turbulence, were not favorable for the creation of the Kepler-223 resonance chain. These initial conditions are indicated using a red font. One can see that, when turbulence is considered, the probability of the successful formation of the correct chain does not differ significantly form all other initial conditions considered, even though they allowed for a correct chain in the nonturbulent case. This suggests that turbulence, in this case, is necessary, while even small turbulence strengths are sufficient to bring the probability in line with other similar initial conditions. The reason for this, however, remains to be investigated in future work.
The trend that less simulations lead to the formation of the Kepler-223 system if the turbulence strength is higher is not apparent for lower values for α, however it can be clearly seen for higher α. One would expect that for an increase in turbulence strength, the amount of simulations that successfully recreate the Kepler-223 resonance chain drops. While this expectation is fulfilled for larger α, it is not as clear for lower α, where simulations with lower surface densities (i.e. lower turbulence strength) seems to not reproduce Kepler-223 as often as those with a higher surface density do. One can see that, once the strength of turbulence lies beyond a critical value, the disruption of the resonance chain is significantly more efficient than below that value. We propose that this is because turbulence below a certain threshold can only disrupt chain formation while the planets are migrating to the correct semi-major axes, while turbulence beyond that threshold can also disrupt chains efficiently after the planets have already fully captured in resonance. This can also serve as an explanation as to why expected trends for the formation of the resonance chain are not as clear for turbulence below that threshold. In the case of only two planets, Batygin & Adams (2017) developed an order-of-magnitude criterion for turbulent disruption of formed resonance chains. It is given by Here, m 1 and m 2 are the masses of the planets, M is the stellar mass, a is the average semi-major axis and f is a factor that accommodates different implementations of the rate of change of the semi-major axis. For our considerations, we apply equation (19) to the second planet pair, because it is observed to be in a 3:2 resonance, which has the lowest resonance index of the chain and should therefore be the least stable. The resulting critical value α crit , obtained by solving (19) toward α, is also indicated in Fig. 3. The analytical values match the results from the simulations in order-of-magnitude. In fact, in the four-planet case, values for α that lie slightly below the predicted α crit from two-planet considerations lead to the disruption of the formed chain. Furthermore, the validity of equation (19) was tested by considering an initial setup where the four planets are already in resonance. Two initial conditions were considered, the first one with a surface density of Σ = 1.39Σ MMSN and = (0.05, 0.02, 0.035), and the second one being Σ = 3.93Σ MMSN and = (0.05, 0.035, 0.035). Equation (19) gives α crit = 6.4 × 10 −2 and α crit = 2.3 × 10 −2 , respectively. Figure 4 shows simulations implementing the aforementioned conditions. It is apparent that the criterion is a good order of magnitude approximation also in the four-planet case, and the behavior of the general case where the planets start outside of the resonances is well-explained by the destruction of formed chains beyond α crit , where the destruction of formed chains becomes possible at lower values for α than α crit because the four-planet case is more sensitive than the two-planet one.

Constraining the minimum turbulence level
While the results from Fig. 3 can give an upper limit for the strength of disk turbulence during the formation of the Kepler-223 resonance chain, it does not provide a lower limit. The amplitude of libration of the Laplace angles can help in finding a lower limit, as Mills et al. (2016) provided the observed amplitudes of these angles.
From an analytical perspective, one would expect the amplitude of the Laplace angles to grow like √ D e ∝ √ αΣ, because the diffusion coefficient governs the average evolution of the eccentricity e, which in turn gives the evolution of the resonance and Laplace angles. The scaling is the same as for the turbulence strength (cfr. equation (13)). For a fixed Σ 0 , this relation can therefore be used to give bounds to α. In Fig. 5, the libration amplitudes for every Laplace angle, grouped by turbulence strength, are shown for the simulations where the correct chain was formed. The figure demonstrates that the libration amplitudes indeed scale like ∼ √ α for constant surface density, allowing the use of the Laplace angle amplitudes to find a lower limit for α. It is shown that for the lowest turbulence strength of α = 10 −4 , the amplitudes of the Laplace angles for the formed chain are smaller than the observed angles. For 10 −3 ≤ α ≤ 10 −2 , an agreement between the observed amplitudes and those resulting from the simulations can be found. For the highest strength of 5 × 10 −2 , the amount of successful simulations is low compared to other strengths, and they are mostly limited to low surface densities. With that in mind, the amplitudes for φ 3 do not agree with the simulations, while the amplitudes for φ 1 and φ 2 can match the expected values. Fig. 3: Amount of successful creations of the Kepler-223 resonance chain for various initial conditions, with varying initial separations and surface densities. For each initial setup, four turbulence strengths α were considered. The green bar denotes successful formations. Generally, a larger turbulence strength leads to fewer successful runs. For each setup, there is a critical value α crit (cfr. equation (19)), for which stronger turbulence significantly hinders the creation of the resonance chain (Batygin & Adams 2017), which is depicted in the plot. This critical value decreases with increasing gas surface density, which is also the case for the amount of successful runs for increasing surface density at constant α (cfr. Eq. 3). The setup marked in a red font did not reproduce the Kepler-223 resonance chain without turbulence, but it does not differ significantly from those that did when the effects of the turbulent disk are considered, indicating that turbulence is a necessary ingredient for this particular setup, which remains to be investigated in future work.

Eccentricities of the planets
Since Mills et al. (2016) also provide the eccentricities of the planets in the Kepler-223 system, these can be compared to the eccentricities that result from the simulations. The comparison between averaged final eccentricity for cases where the correct resonance chain was formed to the eccentricities measured by Mills et al. (2016) is depicted in Fig. 6. The different turbulence strength values that correspond to the successful simulations are indicated by different colors. In the two-planet case, one would not expect the final eccentricities to depend on the turbulence strength. The reason for this is that random kicks as produced by turbulence only lead to a change around the analytical equilibrium point, the value itself is however independent of α and Σ 0 . The analytical value for the two-planet case can be found as (Goldreich & Schlichting 2014) e eq,1 ≈ 1 2(k + 1) where the subscript 1 denotes an evaluation for the inner planet and K is the K-factor as defined in equation (7). This result holds true in the limit of a mass-less inner planet moving outward toward a massive planet on a fixed, circular orbit (circular restricted three-body problem). A more general description shows  Fig. 4: Two setups, also represented in Fig. 3, where turbulence was activated in the simulations only after the formation of the resonance chain. For case (a), equation (19) gives α crit = 6.4 × 10 −2 and for case (b), α crit = 2.3 × 10 −2 . The resulting amount of successful simulations as a function of α is shown, which was compared to the expected critical α crit value (cfr. Eq. 19). The shown results match the theoretical expectations order of magnitude. The positions of the peaks scale with ∼ √ α, as expected due to D e ∝ α, which is the same scaling as for the turbulence strength given by equation (13) for constant surface density Σ 0 . The expected amplitudes according to Mills et al. (2016) are depicted with back dashed lines. The uncertainties of the amplitudes of the Laplace angles are small and therefore are not depicted for our rough value given by those lines. It is apparent that for the lowest turbulence strength α = 10 −4 , the Laplace angle amplitudes are smaller than expected, while for 10 −3 ≤ α ≤ 10 −2 , an agreement between the observed amplitudes and the expected ones can be found. For the highest α = 5 × 10 −2 , the amplitudes resulting from the simulations do not match those expected in the case of φ 3 , while for φ 1 and φ 2 an agreement can still be found.
that a similar relation also holds for the outer planet's equilibrium eccentricity, e eq,2 ∼ K 1/2 2 (Pichierri et al. 2018). Equations (20) and (7) show that the equilibrium eccentricity is proportional to h, and is a function of k. The eccentricity distribution shown in Fig. 6 exhibits multiple peaks, which means that the four-planet case has a dependence on the turbulence strength and initial separation, unlike in the case of only two planets. When comparing our result for the eccentricities with the measured values, we find that there is a mismatch for the two inner-most planets. This could be due to restrictions in our setup, where the inner edge of the protoplanetary disk is modeled by applying high opposing accelerations to the planets close to it, which may not represent the physical situation to high accuracy. In fact, the planet trap could be caused by a sharp drop in surface density, which also affects the eccentricity damping timescale and therefore the equilibrium eccentricity of the planets that are pushed into the trap. By just applying a factor to the acceleration of the semi-major axis, this effect is not considered, resulting in final eccentricities that are lower than they should be. We also note that the eccentricities depend on the aspect ratio, which was fixed at a constant value of h = 0.05. Increasing the aspect ratio of the disk can therefore in principle lead to a better agreement in the eccentricities.

Discussion
In Sect. 3.2.1, we showed that a lower limit for α can be inferred by considering that a certain level of turbulence is needed to reach an excitation of the Laplace angle amplitudes that match those inferred by Mills et al. (2016). While the α values that allow a formation of the Kepler-223 system with the correct Laplace angle libration amplitudes also depend on the disk surface density, a lower limit is given by α ∼ 10 −3 . This value is consistent with turbulence values inferred for an MRI "deadzone" (Flock et al. 2019). Even though α has a space dependence across the protoplanetary disk, the assumption of a constant α is not too strong for our approach, considering that the planetessimals that formed the initial embryos resulting in the Kepler-223 planets formed close to each other, as suggested by embryo growth simulations (e.g., Voelkel et al. 2020Voelkel et al. , 2021. To apply the concept of inferring disk parameters from a planetary system after the disk itself has already dissipated, we adopted the prescription by Cresswell & Nelson (2008) for the additional dissipative forces added to the N-body integrator. This description of planetary migration is comparatively simple, therefore allowing to get gross estimates of parameter ranges and giving a good overview over the idea of this method. For a more thorough investigation, more sophisticated prescriptions for planetary migration in protoplanetary disk may be implemented. While this can of course be achieved by integrating hydrodynamical models that are commonly used for investigations of this kind (see e.g., Bitsch & Kley 2010;McNally et al. 2019;Ataiee & Kley 2020), the goal is to do so with the less computation-intensive method of N-body integration including dissipative forces. In this case, adopting for example models from Paardekooper et al. (2011);Jiménez & Masset (2017), a main aspect is the fact that disk torques responsible for migration in fact exhibit a viscosity dependence, which was not considered here. In addition, such models can include more physically accurate models for the inner edge of the protoplanetary disk, where the positive torque acting as a planetary trap is caused by the slope of the density profile, which changes rapidly close to the inner edge. Adopting these considerations in a self-consistent way including stochastic forces as they are caused by the disk turbulence would allow for a more detailed picture of resonant systems, while in principle applying the same arguments as with the simplified approach.
While an N-body approach is of course not as detailed as a full hydrodynamical consideration, it serves as a good way to reduce the large parameter space, because the N-body approach saves a significant amount of computation time compared to a full hydrodynamical approach. In fact, 6000 simulations were run for 10 5 years of physical time to produce Fig. 3, which corresponds to about 5 × 10 6 orbits of the innermost planet Kepler-223 b. One simulation used, on average, 24 CPU hours, so one orbit required ∼4.8 × 10 −6 CPU hours (∼17 ms). Under the assumption that a typical hydrodynamical approach takes 4 minutes to integrate one orbit, the methods applied here would have required about 3.3 × 10 6 CPU hours (∼38 years) per simulation. The amount of simulations and high integration time that we employed for our considerations would therefore be unfeasible for a hydrodynamical approach.
While planetary systems with long chains of mean-motion resonances are not common due to reasons stated above, Kepler-223 is certainly not the only system with such a peculiar configuration. Another prominent example is the Trappist-1 system, which consists of seven resonant planets. Similarly to the approach taken for Kepler-223, one can infer parameters of the protoplanetary disk in the phase where planetary migration was relevant, even though the disk has already dissipated, by imposing that the disk must allow such a system to form and stay stable during the timescales we observe. This can supplement studies concerning the properties of observable protoplanetary disks, for instance with ALMA (Rafikov 2017;Dullemond et al. 2018;Flaherty et al. 2017Flaherty et al. , 2018. More broadly, any system with a rare orbital configuration, which in turn offers constraints to a protoplanetary disk that would have formed such a rare system, can be used in this way to infer parameter estimates that offer insights on the formation of protoplanetary systems in general.

Conclusions
1. The observed orbital configuration of the Kepler-223 system can be used to infer information about the protoplanetary disk out of which the system formed. This can be achieved without the need of computationally expensive codes that consider the full hydrodynamical picture, but rather by simply employing an N-body integrator with additional dissipative forces mimicking the relevant planet-disk interactions. However, the approach demonstrated in this work should be seen as a proof-of-concept rather than providing detailed results about the disk parameters that governed the late evolution of the Kepler-223 system. This is because this simple model has a number of limits: No evolution of the protoplanetary disk over time is considered, neither is the space-(or time-) dependence of α. Considering these can affect the formation of the chain (Ataiee & Kley 2020) and also the longterm stability (e.g., Pichierri et al. 2018;Pichierri & Morbidelli 2020). Also, one has to acknowledge that this method probes only the evolutionary phase where migration of planets and disk turbulence are the dominating effects. Nevertheless, our method clearly shows that the dynamical constraint from planetary systems in resonance can be used to constrain the disk parameters during planet formation. 2. For a fixed surface density value Σ 0 , this approach can be used to find constraints for the turbulence parameter α, by requiring the formation of the resonance chain under these  (Mills et al. 2016), while the red dashed line denotes the 1σ error. Note that in order to provide better readability, the range of eccentricities depicted on the abscissa is different for the individual plots. Unlike for the two-planet case, a dependence on turbulence strength and initial separation can be found for this case of four planets, resulting in multiple peaks in the distribution. The eccentricities of the two inner most planets do not match the measured values, while for the outer two planets, an agreement can be found. This is probably related to our implementation of the planet trap at the disk's inner edge, which does not capture the entire physical scenario (see Sect. 4). turbulent conditions. This is done under the simplifying assumption that the planets under investigation can be considered to be close to the resonances that they are observed in today. This can be seen as not too restrictive, considering that it is safe to assume that the planets reached an orbital state close to where they are observed at some point in time.
Conversely, under the assumption of a fixed α, the required range of surface densities can be inferred, implying whether the planets of Kepler-223 have formed during the early time of high surface density or the late time of lower surface density. 3. For lower values for α, it is not straightforward to conclude whether the formation of all aspects of the Kepler-223 orbital architecture can be ensured. This is because for weak turbulence, meaning such that captured resonance chains cannot be disrupted, simple analytical insights cannot be put into use and the impact of the turbulent environment is relevant mostly during the formation stage of the chain. The points at which we recover analytical considerations is reached when also an established chain can be disrupted. The limiting parameter α crit can be estimated from two-planet results following Batygin & Adams (2017), as described in equation (19). 4. The values of α that were considered in this work are in the range of α = 10 −4 to 4 × 10 −2 , which is the range found by studies of disk observable with ALMA (Rafikov 2017;Dullemond et al. 2018;Flaherty et al. 2017Flaherty et al. , 2018. These observations are also consistent with hydrodynamical and MRI simulations (e.g., Turner et al. 2014). When just considering the assembly of the correct chain of mean-motion resonances, it is not possible to give a lower limit for α using the methodology employed here. However, since also the amplitudes of the Laplace angles are known for the Kepler-223 system (Mills et al. 2016), an additional constraint can be provided. As shown in Fig. 5, a certain level of turbulence is needed to reproduce the correct amplitudes of the Laplace angles, which is about 15°for φ 1 and 30°for φ 2 and φ 3 , peak-to-peak. As a lower limit for α, the consideration of the Laplace angle amplitudes gives α ∼ 10 −3 , which is consistent with values typically considered for an MRI "dead-zone" (Flock et al. 2019). Our measurement of viscosity comes with the uncertainty of an unknown value for the surface density, since our simulations allow us to make statements about Σ 0 √ α. For example, a higher turbulent viscosity can be consistent with our results if we adopt a lower surface density. It was also found that even if the correct MMR chain is formed, the system can still fail to reproduce the correct amplitudes for the Laplace angles and can therefore not be considered to have correctly reproduced the observations, mainly important at high α. 5. While the described methods can be used to find rough constraints for disk parameters during the later stages of planet formation and the assembly of planetary systems, it has shortcomings, which can be seen by the example of the final planet eccentricities as depicted in Fig. 6. The eccentricities of the inner two planets as measured by Mills et al. (2016) were not reproduced here. However, this does not invalidate our method, because on the one hand, a global mismatch of the eccentricities can be fixed by modifying the aspect ratio prescription, where an increase could lead to a shift in eccentricities that provides a better match with the observed values (Goldreich & Schlichting 2014;Pichierri et al. 2018).
On the other hand, we propose as a reason for the presented mismatch the model of the inner edge of the protoplanetary disk, which is strongly simplified and does not accurately represent all physical properties. While the focus of this investigation was on the methodology itself, this issue could be fixed by employing more physically accurate models (e.g., Izidoro et al. 2017Izidoro et al. , 2019).
Our method clearly shows that the currently observed system architectures hold imprints to their formation history. Especially systems with multiple planets in resonance can help to constrain disk parameters. This is especially important to also probe different environments, like in the Trappist-1 system, orbiting an Mdwarf, where constraints from disk observations are even rarer.