Martin Jutzi
email: martin.jutzi@space.unibe.ch

Patrick Michel

Collisional heating and compaction of small bodies: Constraints for their origin and evolution

Keywords: Asteroids, collisions, Collisional physics

published or not. The documents may come

Collisional heating and compaction of small bodies:

Constraints for their origin and evolution

Introduction

The various objects composing the different small body populations, such as asteroids and comets, are the results of a billion-year-long collisional evolution of varying intensity. Their current properties provide important clues to their origin and history, and therefore represent a window to the early stages of the formation of the Solar System. Most asteroids smaller than about 50 km in diameter are the result of a break-up of a larger parent body [START_REF] Bottke | The fossilized size distribution of the main asteroid belt[END_REF]. They have been modified by subsequent collisions and their shapes, interior structure and spin state are determined to a large degree by the last major (global scale) impact event. Collisional processes -ranging from low-velocity mergers to catastrophic disruptions -have been proposed to also play an important role in the evolution of Kuiper Belt Objects and cometary nuclei (e.g. [START_REF] Durda | Collision Rates in the Present-Day Kuiper Belt and Centaur Regions: Applications to Surface Activation and Modification on Comets, Kuiper Belt Objects, Centaurs, and Pluto-Charon[END_REF][START_REF] Morbidelli | Comets as collisional fragments of a primordial planetesimal disk[END_REF], although they are generally less frequent and less intensive than collisions between asteroids.

Important questions related to collisions are: at which level were small bodies processed by past collisions and to what extent do they retain a record of processes that took place during the formation and early evolution of the Solar System? By analysing the degree of material processing (e.g. heating and compaction) during impacts, constraints for the origin and evolution of the small body populations as a whole as well as of individual objects -such as the recently visited asteroids Ryugu and Bennu or comet 67P/ChuryumovGerasimenko -can be determined. Here, we perform such an analysis from a large data set of small body collisions produced in a previous study in which the focus was on the size and ejection speed distributions of resulting fragments [START_REF] Jutzi | Fragment properties from large-scale asteroid collisions: I: Results from SPH/N-body simulations using porous parent bodies and improved material models[END_REF], and not on thermodynamical and physical state aspects.

Previous studies of impact heating in small body collisions (e.g. [START_REF] Davison | Numerical modelling of heating in porous planetesimal collisions[END_REF][START_REF] Davison | The early impact histories of meteorite parent bodies[END_REF][START_REF] Ciesla | Thermal consequences of impacts in the early solar system[END_REF] have considered a limited range of conditions, mainly focusing on head-on impacts. The effect of impact obliquity was studied by [START_REF] Davison | The effect of impact obliquity on shock heating in planetesimal collisions[END_REF]. However, these studies did not distinguish between heated material bound to the largest remnant and the heated material which is escaping the system. Recently, [START_REF] Wakita | Fates of hydrous materials during planetesimal collisions[END_REF] performed head-on impacts onto non-porous objects to specifically investigate the fate of hydrous materials. The simulations of small body break-ups by [START_REF] Jutzi | Fragment properties from large-scale asteroid collisions: I: Results from SPH/N-body simulations using porous parent bodies and improved material models[END_REF], which are analysed here, used porous targets of 50 -400 km in diameter and investigated a large range of impact velocities, angles and energies. We investigate the heating and compaction produced in these impacts and study how much processed material is retained (i.e., is bound to the largest remnant).

In a recent study by [START_REF] Schwartz | Catastrophic disruptions as the origin of bilobate comets[END_REF], it was found that not much heating or compaction occurs during the disruption of comets of 10 km in size, except for the small ejecta that do not contribute to the reaccumulation process producing the hierarchy of aggregates composing the main fragment distribution. Analysing the degree of material processing for bodies of 50 -400 km in size allows us to understand how such a processing depends on scale and to apply it to small comet formation as well as to the formation of asteroid families whose parent bodies were in this size range, with implications for their physical properties. For instance, the NASA mission OSIRIS-REx and the JAXA mission Hayabusa2 both found that their asteroid targets, Bennu and Ryugu, have an hydration band, which is much deeper for Bennu than for Ryugu [START_REF] Hamilton | Evidence for widespread hydrated minerals on asteroid (101955) Bennu[END_REF][START_REF] Kitazato | The surface composition of asteroid 162173 Ryugu from Hayabusa2 near-infrared spectroscopy[END_REF]. It is thought that they come from asteroid families whose parent body is in the 100 km size range [START_REF] Walsh | Introducing the Eulalia and new Polana asteroid families: Reassessing primitive asteroid families in the inner Main Belt[END_REF]. When such an asteroid family is created, does the material experience enough heating that any original hydrous material would be lost? Answering this question puts important constraints on the origin and evolution of small hydrated asteroids. In the context of comet formation, the analysis of heating and compaction in break-ups of 50 -400 km objects allows us to address the question whether or not comets were born big [START_REF] Youdin | Streaming Instabilities in Protoplanetary Disks[END_REF][START_REF] Johansen | Rapid planetesimal formation in turbulent circumstellar disks[END_REF][START_REF] Cuzzi | Towards initial mass functions for asteroids and Kuiper Belt Objects[END_REF][START_REF] Morbidelli | Asteroids were born big[END_REF], and what is their maximal initial size. Can small comets such as 67P originate from the disruption of such large bodies, or must they come from smaller ones that experience less heating as found in previous studies (Jutzi and Benz, 2017;[START_REF] Schwartz | Catastrophic disruptions as the origin of bilobate comets[END_REF])?

Here we concentrate on the processing of 1) material that remains bound to the largest remnant of a collisional disruption and 2) all the escaping material. In a future study, we will analyze the formation of aggregates and the fraction of heated and compacted material in these aggregates using a N-body code as in [START_REF] Schwartz | Catastrophic disruptions as the origin of bilobate comets[END_REF]. This will allow for a more complete understanding of how heated and compacted materials are distributed in the whole set of fragments (reaccumulated and escaping ones) produced by a large disruption. The possible modification of the interiors by early thermal evolution is investigated in a separate study (Golabek and Jutzi, in preparation).

In Sec. 2, we outline our modeling approach and in Sec. 3, we present the results of our analyses of the whole set of impact simulations covering a wide range of impact energies and angles. The implications of the results are discussed in Section 4.

Model approach

The numerical scheme that was used to produce the data set of simulations that is analysed here is described in detail in paper I [START_REF] Jutzi | Fragment properties from large-scale asteroid collisions: I: Results from SPH/N-body simulations using porous parent bodies and improved material models[END_REF].

Here we provide a short summary of our collision modeling and present our approach to calculate the degree of heating and compaction in the largest collisional remnant and the escaping fragments.

Collision modeling

To model an asteroid disruption and subsequent reaccumulations, [START_REF] Jutzi | Fragment properties from large-scale asteroid collisions: I: Results from SPH/N-body simulations using porous parent bodies and improved material models[END_REF] used an SPH / N-body approach as introduced by [START_REF] Michel | Collisions and Gravitational Reaccumulation: Forming Asteroid Families and Satellites[END_REF][START_REF] Michel | Disruption of fragmented parent bodies as the origin of asteroid families[END_REF]. The early phases of the collisions were simulated using an SPH impact code (Benz andAsphaug, 1994, 1995;[START_REF] Nyffeler | Modelling of Impacts in the Solar System on a Beowulf Cluster[END_REF][START_REF] Jutzi | Numerical simulations of impacts 19 involving porous bodies. I. Implementing sub-resolution porosity in a 3D SPH hydrocode[END_REF][START_REF] Jutzi | SPH calculations of asteroid disruptions: The role of pressure dependent failure models[END_REF] that includes self-gravity as well as material models including strength, friction and porosity. The N-body code pkdgrav was then used to compute the dynamical evolution of the system to late time. [START_REF] Jutzi | Fragment properties from large-scale asteroid collisions: I: Results from SPH/N-body simulations using porous parent bodies and improved material models[END_REF] used porous parent bodies with a nominal diameter of D = 100 km. To study the scale-dependence of the results, they additionally investigated D = 50, 150 and 400 km targets. For the D = 100 km objects, the same matrix of impact conditions as used in [START_REF] Durda | Size-frequency distributions of frag-ments from SPH -N-body simulations of asteroid impacts: Comparison with observed asteroid families[END_REF] and [START_REF] Benavidez | A Comparison Between Rubble-Pile and Monolithic Targets in Impact Simulations: Application to Asteroid Satellites and Family Size Distributions[END_REF] was explored covering a wide range of impact speeds (from 3 to 7 km/s), impact angles (from 15 • to 75 • with 15 • increments)1 and impactor diameters, allowing a comparison with those two studies. In these studies, reaccumulating N-body particles were simply merged, as the main focus was on the fragment size and ejection speed distributions. The resulting shapes, the compaction of material and the distribution of heat were not investigated.

Here we focus on the degree of heating and compaction produced by the impact, which are all part of the outcome of the SPH simulations. In our analysis, we distinguish between the largest remnant and material not bound to this remnant. The properties of the individual smaller fragments will be investigated in a future study.

Determination of largest remnant and unbound fragments

We analyse the simulation data obtained by [START_REF] Jutzi | Fragment properties from large-scale asteroid collisions: I: Results from SPH/N-body simulations using porous parent bodies and improved material models[END_REF] in a postprocessing step in terms of heating and compaction of the largest remnant and the escaping fragments. The simulations performed in [START_REF] Jutzi | Fragment properties from large-scale asteroid collisions: I: Results from SPH/N-body simulations using porous parent bodies and improved material models[END_REF] did not track the particle IDs in the N-body calculations, which would be required to follow the individual history of the fragments. Therefore, we compute the mass of the largest remnant M lr and the corresponding mass of the unbound fragments M f at the end of the impact phase using an iterative energy balance approach [START_REF] Benz | Catastrophic Disruptions Revisited[END_REF][START_REF] Jutzi | Fragment properties at the catastrophic disruption threshold: The effect of the parent body's internal structure[END_REF]. This approach was demonstrated in previous studies to be consistent with the outcome of N-body simulations of the gravitational phase, as long as M lr is larger than a few 10% [START_REF] Jutzi | Fragment properties at the catastrophic disruption threshold: The effect of the parent body's internal structure[END_REF] of the total mass. We apply this scheme to the simulations with specific impact energies Q smaller than 2 times the catastrophic disruption threshold Q * D . For the small subset of simulations which have higher energies, we do not distinguish between largest remnant and smaller fragments.

Heating analysis

The shock physics code calculations of the collisions were performed in [START_REF] Jutzi | Fragment properties from large-scale asteroid collisions: I: Results from SPH/N-body simulations using porous parent bodies and improved material models[END_REF] using the relatively simple Tillotson equation of state which does not allow for a direct computation of a thermodynamically consistent temperature. However, we can still approximately compute the temperature increase dT caused by an impact by analysing the increase in specific internal energy. As shown recently, this approach leads to temperatures comparable to the ones computed via the more sophisticated ANEOS equation of state, as long as vaporisation is negligible [START_REF] Emsenhuber | SPH calculations of Marsscale collisions: The role of the equation of state, material rheologies, and numerical effects[END_REF]. In [START_REF] Schwartz | Catastrophic disruptions as the origin of bilobate comets[END_REF], it was assumed that the heat capacity is constant, while in reality it depends on temperature. To improve the realism, here we use a temperature dependent heat capacity to compute the temperature increase for rocky as well as cometary-like materials:

T = T 0 + ∆u i /c(T ) (1)
where T 0 is the initial temperature, ∆u i the increase in specific energy and c(T ) the heat capacity. We then take into account the temperature depen-dence of the heat capacity by using the simple relation:

c(T ) = b + aT (2)
defined by the parameters a and b. Using the equations ( 1) and ( 2) we obtain the following relation for the temperature as a function of specific internal energy:

T = T 0 + ∆u i c(∆u i ) (3) with c(∆u i ) = 1 2 T 0 a + 1 2 b + 1 2 T 2 0 a 2 + 2T 0 ab + b 2 + 4∆u i a (4)
For rocky asteroids, we compute c rock = c(∆u i ) with a = 2.5 and b = 0, which provides a rough approximation to laboratory measurements of the heat capacity of forsterite at temperatures below 400 K [START_REF] Robie | Heat capacities and entropies of Mg2SiO4,Mn2SiO4,and Co2SiO4 between 5 and 380 K[END_REF].

These parameters, together with T 0 = 50 K, result in c rock = 1000 J/(kg K) at T = 400 K (corresponding to ∆u i = 3.5 10 5 J/kg). For T > 400 K, we use a constant c rock = 1000 J/(kg K).

For our analysis of cometary-like material we use a weight-averaged heat capacity [START_REF] Davidsson | The primordial nucleus of comet 67P/Churyumov-Gerasimenko[END_REF]:

c mix = f rock c rock + f ice c ice (5)
where

f ice = 1 -f rock , f rock = m r,rock-ice /(1 + m r,rock-ice
) and we assume a rock-to-ice mass ratio of m r,rock-ice = 2. The heat capacity of water ice c ice is computed using the parameters a = 7.49, b = 90 Ws/kg/K [START_REF] Klinger | Some consequences of a phase transition of water ice on the heat balance of comet nuclei[END_REF].

In this analysis we also compute the mass fraction of sublimated water-ice for T > 180 K as

f = (∆u i -u crit )/u vap (6)
where u crit = 188 kJ/kg is the specific energy at T = 180 K (using T 0 = 50 K) and u vap = 3000 kJ/kg is the heat of vaporization [START_REF] Ahrens | Shock Vaporization and the Accretion of the Icy Satellites of Jupiter and Saturn[END_REF].

Compaction

In the collision calculations by [START_REF] Jutzi | Fragment properties from large-scale asteroid collisions: I: Results from SPH/N-body simulations using porous parent bodies and improved material models[END_REF], porosity was modeled using a sub-resolution approach based on the P-alpha model [START_REF] Jutzi | Numerical simulations of impacts 19 involving porous bodies. I. Implementing sub-resolution porosity in a 3D SPH hydrocode[END_REF]. The material properties (crush-curve) of the porous target used in this study were those that provided the best match to impact experiments on pumice targets [START_REF] Jutzi | Numerical simulations of impacts involving porous bodies. II. Comparison with laboratory experiments[END_REF]. The porosity model takes into account the enhanced dissipation of energy during compaction of porous materials.

In order to compute the degree of compaction caused by the collision, we consider the relative change of porosity, defined as

p 0 -p p 0 = α 0 -α α(α 0 -1) (7)
where p 0 = 1 -1/α 0 is the initial porosity, p = 1 -1/α is the post-impact porosity, and α 0 and α are the initial and post-impact distention, respectively.

For full compaction, (p 0 -p)/p 0 = 1 while in the case of no compaction, (p 0 -p)/p 0 = 0.

Results

Impact heating of rocky material

The degree of impact heating of rocky material computed according to Eq. ( 3) is shown in Fig. 1 for D = 100 km targets as a function of the normalized specific impact energy for three different impact velocities and two impact angles (15 • and 45 • ). For this computation we use the post-impact state of the largest remnant and the unbound material, 400 s after the impact, i.e. when the shock wave has long passed (see section 3.1.1). The temperature increase is computed via Eq. (3) using c rock . Also displayed by the black and gray curves are the mass fraction of the largest remnant M lr /M tot , the mass fraction of the escaping material M esc /M tot = 1 -M lr /M tot as well as the fraction of thermal energy contained in the largest remnant and the ejected material, respectively.

We find that the fraction of thermal energy in the largest remnant is always (much) smaller than the mass fraction of the largest remnant, while

for the escaping material it is the opposite. In other words, there is a clear dichotomy in terms of impact processing: the escaping material always experiences stronger heating than the material bound to the largest remnant. It is worth pointing out that this is not only the case for small ratios of Q/Q * D (cratering regime), but also in catastrophic or super-catastrophic collisions,

where more than half of the initial mass is dispersed. As can be seen in Fig. 1, the differences are larger for more oblique impact angles.

Interestingly, for a given impact angle and velocity, the temperature increase experienced by the escaping material is roughly constant, as larger projectiles heat a larger region of the target, but at the same time eject a larger amount of target material. The heating of the bound material (largest remnant) increases with increasing impact energy, but is less than the heating of the escaping material. For the bound material, there is a maximum heating (which depends on impact angle and velocity) because the relative fraction of thermal heat contained in the largest remnant decreases faster with impact energy than its relative mass fraction (gray vs. black lines with filled symbols in Fig. 1). In other words, a larger fraction of heated material escapes at high impact energies. We also compute the mass fractions experiencing a specific temperature increase of dT > 400 K for a range of impact parameters (Fig. 2 and 3), which confirms these observations.

Convergence tests

To quantify the dependence of our results on numerical parameters, the degree of impact heating and the mass of largest remnant are shown in Fig. 4 as a function of numerical resolution for the case with 45 • impact angle and velocity of 3 km/s using different projectile sizes and corresponding specific impact energies.

The dependence of the results on the time of the analysis is shown in Fig. 5.

Overall, there is a reasonable convergence of the investigated quantities, especially given the potentially much larger effects of material properties, impact geometries, target structures and shapes, etc.

Compaction

The degree of compaction (Fig. 6) shows very similar characteristics as the degree of heating. The escaping material experiences a relatively high degree of compaction, independent of the impact energy, while the compaction of the bound material increases with increasing impact energy, but is much less significant. Figure 7 shows the mass fraction of material compacted by more than 50% for specific cases. What is not taken into account in this analysis is the addition of macro-porosity by the reaccumulation process, which can again increase the overall porosity (as studied in Jutzi and Benz ( 2017)).

Impact heating of cometary-like material

We also compute the temperature increase that a cometary-like material would experience due to impacts. To do this we use the heat capacity for a rock-ice mixture as defined in Eq. ( 5). Figure 8 shows the mass fractions of material heated by dT > 80 K for the cases with 3 km/s impact velocity and varying impact angles.

Most of the simulations in [START_REF] Jutzi | Fragment properties from large-scale asteroid collisions: I: Results from SPH/N-body simulations using porous parent bodies and improved material models[END_REF] were performed using D = 100 km targets. However, a subset of impact parameters were studied using a range of targets sizes (50 -400 km). These runs are used here to investigate the size dependence of the degree of heating and compaction. The results are summarised in Table 4.2 and displayed in Figs. 9 and 10.

The mass fraction of sublimated water-ice computed according to Eq. ( 6) is shown in Fig. 11. Except in the cases of the most energetic, supercatastrophic collisions, the impacts cause only little vaporization.

Discussion and perspectives

The results of our analysis of compaction and heating effects in collisional break-ups of 50 -400 km sized bodies have important implications regarding the formation (initial size) and evolution of small solar system bodies.

Constraints for the formation of asteroids Ryugu and Bennu

The asteroids Rygu and Bennu have been suggested to have originated from the Polana/Eulalia family forming events [START_REF] Walsh | Introducing the Eulalia and new Polana asteroid families: Reassessing primitive asteroid families in the inner Main Belt[END_REF][START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF]. With parent body sizes in the range of D = 100 km, the masses of these subkm-sized objects are only ∼ 10 -6 of the parent body mass. Current hybrid SPH-N-body simulations can therefore not properly resolve explicitly the formation of such small fragments (reaccumulated rubble piles), as they would only be represented by only 1-10 SPH particles in the original D ∼ 100 km parent body. However, our results (Fig. 2 and3)

suggest that the different degree of hydration observed on Ruygu and Bennu [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF][START_REF] Lauretta | The unexpected surface of asteroid (101955) Bennu[END_REF] does not necessarily mean that they come from different parent bodies, or that they experienced different degrees of dehydration after their formation, due for instance to their different perihelion histories. It could rather be due to the different degree of heating of the material that reaccumulated to form both of them in a single disruption event. In effect, depending on impact velocity and angle, up to 20% of unbound material is heated > 400 K (i.e., Ryugu could have been reaccumulated from this material), while the rest is heated less (i.e.

this could be Bennu forming material). We note that given the small size of Ryugu and Bennu compared to the parent body, they are likely sourced from localized regions, which would be consistent with the uniform degree of hydration suggested by the spectral observations.

A more detailed analysis investigating the degree of heating in individual fragments and a discussion in the context of the formation of Ryugu and Bennu, including their shape, is presented in [START_REF] Michel | Formation of top-shaped asteroids by collisional disruption: implications for the origin of Ryugu and Bennu[END_REF].

Contraints for the formation of comeatry nuclei

Previous studies (e.g. Jutzi and Benz, 2017;Jutzi et al., 2017;[START_REF] Schwartz | Catastrophic disruptions as the origin of bilobate comets[END_REF] found that at small < 10 km scales, even catastrophic disruptions do not cause significant heating of the largest fragments. The collisions be-tween ∼ 100 km sized objects analysed here take place at significantly larger scales. Because of the size-dependence of the catastrophic disruption energy, there is significantly increased heating (as Q * D strongly increases with increasing size). According to the scaling law Q * D = aR 3µ V2-3µ using µ = 0.42 (Jutzi et al., 2017) and a constant velocity, Q * D,R=50km /Q * D,R=3.5km ∼ 30, we can expect an overall increase of heating (compared to [START_REF] Schwartz | Catastrophic disruptions as the origin of bilobate comets[END_REF] by the same factor. Indeed, our results (e.g., Fig. 8) show significantly more heating than found by [START_REF] Schwartz | Catastrophic disruptions as the origin of bilobate comets[END_REF] at smaller scales for comparable impact velocities (up to 1 km/s) and impact angles 2 . However, as discussed above, there is also a clear dichotomy in the degree of heating between the bound and unbound material.

For instance in the case of a 3 km/s impact with 45 • angle and Q/Q * D ∼ 1.1, we find that ∼ 4% of the bound material experiences a temperature increase dT > 80 K (i.e., 96% is heated by less than 80 K). On the other hand, around 60% of the unbound material experiences a temperature increase dT > 80 K. For less disruptive impacts in the cratering regime with Q/Q * D ∼ 0.01, we find that only ∼ 1% of the bound material experiences a temperature increase dT > 80 K, while still ∼ 40% of unbound material has dT > 80 K.

The analysis of the size dependence of these processes in collisional disruptions (Fig. 9 and 10 ; Table 4.2) indicates that for bodies smaller than about R = 25 km, the heating as well as the compaction of the bound material become negligible. This may suggest that the maximal initial size of cometary nuclei -assuming that comets were born big -must have been in this size range, in order for small comets resulting from the disruption of these bodies to keep their primitive properties.

However, the distribution of heat as a function of fragment size will require a more detailed analysis, following explicitly the evolution of the unbound material with a N-body code, as done by [START_REF] Schwartz | Catastrophic disruptions as the origin of bilobate comets[END_REF] and [START_REF] Michel | Formation of top-shaped asteroids by collisional disruption: implications for the origin of Ryugu and Bennu[END_REF], which will be performed more generally in a subsequent study. Another important aspect is the possible modification of the comet precursors' interiors by early thermal evolution (Golabek and Jutzi, in preparation).

The complete understanding of small body heating needs to account for many aspects, and our study on impact heating and compaction, indicating in which case such processes are expected, provides direction for future investigations. and compacted by > 50% as a function of normalized impact energy Q/Q * D . The data corresponds to the simulations with a 45 • impact angle, a velocity of 3 km/s and varying target radii R.

R (km) State Q/Q * D f >80K f >50% 25 bound 1.08E-

Figure 1 :

 1 Figure 1: Degree of heating in simulations using D = 100 km targets, various impact angles (15 • and 45 • ) and velocities (3, 5, 7 km/s). In each case the mass fraction of material heated above a certain temperature for bound material (left column), all material (middle column) and escaping fragments (right column) is shown. The left and right columns are slightly shifted for a better visualisation. The x-axis corresponds to the specific impact energy Q normalised by the catastrophic disruption threshold Q * D (which is computed based on the data by Jutzi et al. (2019) for each combination of impact angle and velocity). For Q/Q * D > 2, only the middle column is shown. The black line with filled symbols indicates the mass fraction of the largest remnant, while the black line with open symbols corresponds to the mass fraction of the escaping material. The gray lines show the corresponding fraction of thermal heat (filled symbols: contained in largest remnant; open symbols: contained in escaping material).

Figure 2 :

 2 Figure 2: Fraction of material heated by dT > 400 K. Shown are the simulations with 45 • impact angle and velocities of 3, 5 and 7 km/s. The gray lines show the relative sizes of the largest remnants.

Figure 3 :

 3 Figure 3: Fraction of material heated by dT > 400 K. Shown are the simulations with varying impact angles and a velocity 5 km/s. The gray lines show the relative sizes of the largest remnants.

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Fraction of unbound material heated by dT > 400 K and relative size of largest remnant (M lr ). Shown are the simulations with 45 • impact angle and velocity of 3 km/s as a function of impact energy, using different numerical resolutions: 100 k, 400 k (nominal) and 1000 k SPH particles.

Figure 7 :

 7 Figure 7: Fraction of material compacted by more than 50%. Shown are simulations with various impact angles and a velocity of 3 km/s.

Figure 8 :

 8 Figure 8: Fraction of material heated by dT > 80 K using cometary-like properties. Shown are simulations with various impact angles and a velocity of 3 km/s.

Figure 9 :

 9 Figure9: Fraction of material heated by dT > 80 K using cometary-like properties. Shown are the simulations with a 45 • impact angle, a velocity of 3 km/s and varying target sizes.

Table 1 :

 1 Fraction f of material heated by dT > 80 K (using cometary-like properties)

			02 4.43E-03 4.81E-03
	25	bound	2.76E+00 1.00E-04 5.76E-03
	25	unbound 1.08E-02 2.50E-01 2.29E-01
	25	unbound 2.76E+00 4.19E-01 5.83E-01
	50	bound	4.21E-03 5.11E-03 5.59E-03
	50	bound	1.09E-02 1.12E-02 1.52E-02
	50	bound	2.77E-02 2.12E-02 3.32E-02
	50	bound	6.95E-02 3.53E-02 6.19E-02
	50	bound	1.74E-01 4.75E-02 9.34E-02
	50	bound	4.41E-01 4.82E-02 1.14E-01
	50	bound	1.11E+00 3.83E-02 1.20E-01
	50	unbound 4.21E-03 4.23E-01 3.34E-01
	50	unbound 1.09E-02 4.32E-01 4.65E-01
	50	unbound 2.77E-02 5.07E-01 5.82E-01
	50	unbound 6.95E-02 5.43E-01 6.52E-01
	50	unbound 1.74E-01 5.68E-01 7.04E-01
	50	unbound 4.41E-01 5.79E-01 7.40E-01
	50	unbound 1.11E+00 6.10E-01 7.88E-01
	100	bound	1.74E-03 5.73E-03 6.34E-03
	100	bound	4.43E-01 1.46E-01 2.96E-01
	100	unbound 1.74E-03 6.79E-01 4.76E-01
	100	unbound 4.43E-01 8.09E-01 9.27E-01
	200	bound	6.95E-04 5.70E-03 3.38E-01
	200	bound	1.77E-01 2.25E-01 6.81E-01
	200	unbound 6.95E-04 9.88E-01 7.33E-01
	200	unbound 1.77E-01 8.84E-01 9.82E-01

We note that for a random impact in 3 dimensions, the probability function dP for the impact angle θ is given by dP = sin(2θ)dθ, giving a maximum probability of impact at θ = 45 • . For a discussion of quasi

2-dimensional systems (e.g. flat disks), see[START_REF] Leleu | The peculiar shapes of Saturn's small inner moons as evidence of mergers of similar-sized moonlets[END_REF]

We note that[START_REF] Schwartz | Catastrophic disruptions as the origin of bilobate comets[END_REF] used the same impact code, but slightly different material parameters, such as a lower density and crushing pressures.

Acknowledgments

The authors acknowledge funding support from the European Union's Horizon 2020 research and innovation programme under grant agreement No 870377 (project NEO-MAPP). M.J. acknowledges support from the Swiss National Centre of Competence in Research PlanetS. P.M. acknowledges funding support from the French space agency CNES as well as from Academies of Excellence: Complex systems and Space, environment, risk, and resilience, part of the IDEX JEDI of the Université Côte dAzur. We also thank two anonymous reviewers for their constructive and valuable comments.


S Abadian
email: sevag.abadian@c2n.upsaclay.fr

G Magno

V Yam

B Dagens

INTEGRATED MAGNETO-PLASMONIC ISOLATION ENHANCEMENT BASED ON COUPLED RESONANCES IN SUBWAVELENGTH GOLD GRATING

Keywords: magnetoplasmonic, isolators, waveguides, gold grating, nanoslits, Fabry Perot resonance, telecom wavelength

In this paper, we propose and numerically investigate a new non-reciprocal optical device enabling enhanced optical TE mode isolation at telecom wavelength, compatible with applications. The principle is based on magnetoplasmonic interaction between a garnet waveguide and a side-coupled gold grating. Grating nanoslits host Fabry-Perot resonances that couple to Surface Plasmon Polaritons excited at the interface between the metal and the waveguide. Additionally, coarse or fine sensitivity of the Fabry Perot resonances to nanoslits dimensions opens the way to continuously selectable operation frequency of the plasmonic isolator.

Introduction

Magnetoplasmonics is a paradigm which benefits from the interaction between surface plasmon polaritons (SPPs) and magneto-optical (MO) materials. On the one hand, plasmonics is employed to attain extreme electromagnetic field enhancement and subwavelength localization. On the other hand, MO materials bring optical non-reciprocity properties [START_REF] Armelles | Magnetoplasmonics: Combining Magnetic and Plasmonic Functionalities[END_REF]. The transverse magnetooptical Kerr effect (TMOKE) can be fruitfully exploited in compact integrated guiding structures to achieve optical isolation when the MO material magnetization is perpendicular to the plane formed by the direction of light propagation and its polarization. In this case, a change in the intensity and phase of the reflected or transmitted light at the MO material interface occurs without modification of its polarization [START_REF] Zvezdin | Modern Magnetooptics and Magnetooptical Materials[END_REF]. The design of optical guided structure with non-reciprocal behavior requires simultaneous spatial and time-reversal symmetry breaking in the waveguide. Spatial symmetry breaking is easily achieved by stacking layers of different refractive indices, whereas time-symmetry breaking can be provided by the magneto-optical material [START_REF] Stadler | Integrated Magneto-Optical Materials and Isolators: A Review[END_REF]. Optical isolators can remove optical feedback into laser diodes which could disturb the active medium, thus making them necessary components in telecommunication systems. More generally the achievement of a large magneto-optic response is essential for many applications like optical modulators [START_REF] Sobolewski | Magneto-optical modulator for superconducting digital output interface[END_REF], switches [START_REF] Bahuguna | Magneto-Optic-Based Fiber Switch for Optical Communications[END_REF], optical isolators [START_REF] Shoji | Silicon Mach-Zehnder interferometer optical isolator having 8 nm bandwidth for over 20 dB isolation[END_REF] or circulators [START_REF] Ghosh | Adhesively bonded Ce:YIG/SOI integrated optical circulator[END_REF], mode-conversion waveguides and optical data storage [START_REF] Hansen | Magneto-optical recording materials and technologies[END_REF].

Many isolator designs have been proposed, along with different materials and different geometries. Vanwolleghem et al. [START_REF] Vanwolleghem | Design and optimization of a monolithically integratable InP-based optical waveguide isolator[END_REF] proposed a TM waveguide isolator with a transversely magnetized ferromagnetic metal on the top. A similar structure based on the TE mode was presented by Shimizu and Nakano [START_REF] Shimizu | Fabrication and Characterization of an InGaAsp/InPActive Waveguide Optical Isolator With 14.7 dB/mm TE Mode Nonreciprocal Attenuation[END_REF] with vertically magnetized MO material on the side of the waveguide. In both cases, TMOKE induces a nonreciprocal shift of the complex effective index of the guided mode resulting in an absorption which is different in forward and backward propagation directions [START_REF] Vanwolleghem | Design and optimization of a monolithically integratable InP-based optical waveguide isolator[END_REF][START_REF] Shimizu | Fabrication and Characterization of an InGaAsp/InPActive Waveguide Optical Isolator With 14.7 dB/mm TE Mode Nonreciprocal Attenuation[END_REF]. However, optical losses which are associated with ferromagnetic metals are high and the idea of substituting them by low-absorptive magneto-optical dielectrics like garnet was introduced. When considered on their own, these lossless MO materials do not lead to any significant optical intensity transmission nonreciprocity due to the TMOKE, and are used for the non-reciprocal phase shift they induce in interferometric or resonant structures [START_REF] Stadler | Integrated Magneto-Optical Materials and Isolators: A Review[END_REF]. To enhance the impact of the TMOKE, these MO materials can be incorporated with plasmonic structures which can induce SPP excitation into the interferometric system [START_REF] Firby | Magnetoplasmonic Isolators Utilizing the Nonreciprocal Phase Shift[END_REF]. Non-reciprocal losses also are induced by the combination of plasmonic losses and magneto-optical effect, and in addition they can be enhanced by plasmonic metal nanostructuration [START_REF] Belotelov | Extraordinary transmission and giant magneto-optical transverse Kerr effect in plasmonic nanostructured films[END_REF][START_REF] Halagacka | Magnetoplasmonic nanograting enables optical nonreciprocity sign reversal[END_REF]. Nevertheless, these properties have never been used in integrated isolators.

We consider here a new isolator configuration consisting of a MO waveguide side coupled to a nanostructured gold grating which could be efficiently integrated in Photonic Integrated Circuits (PICs). Here the grating period is subwavelength in order to avoid Bragg reflections in the waveguide which could deteriorate the isolation performance. Such a system benefits from the intertwined application of the TMOKE, since this effect is well-adapted to the waveguiding configuration which is based on total internal reflection on opposing interfaces, and the joint plasmonic resonances inside the gold grating. This arrangement is ensured to break spatial and time symmetry simultaneously and permits coupling between the SPPs, the waveguide mode and the horizontal and vertical Fabry-Perot resonant modes in the grating slits [START_REF] Jiao | Fabry-Perot like phenomenon in the surface plasmons resonant transmission of metallic gratings with very narrow slits[END_REF]. Thus, the TMOKE, which when used apart is a weak effect, can gain interest and act a key role in developing new PIC components when adequately enhanced by SPPs. Finally, the isolation functionality of the proposed structure depends on the coupling degree of the fundamental TE mode in the MO waveguide to the SPPs propagating at the metal-dielectric interface and the excited metal-insulator-metal (MIM) modes inside the slits. As we show in this paper, this mechanism contributes to the enhancement of non-reciprocal transmission through the MO waveguide bringing us closer to finding integrable optical signal isolators.

Structure and principle

The proposed structure, illustrated in Fig. 1, mainly consists of a magneto-optical waveguide and a noble metal grating. The MO Bismuth Iron Garnet (𝐵𝐼𝐺) ridge waveguide has dimensions of height ℎ 𝐵𝐼𝐺 = 0.45 𝜇𝑚 (y direction) and width 𝑤 𝐵𝐼𝐺 = 1 𝜇𝑚 (z direction) chosen in order to allow the propagation of the fundamental guided Transverse Electric (TE) mode. The confinement factor 𝛤, defined as the proportion of the power of a given mode that lies within the core waveguide, is around 60% at wavelength 1.55 𝜇𝑚.

The waveguide is embedded in a dielectric medium with a dielectric constant 𝜀 𝑆𝑖𝑂 2 = 1.45 2 and placed on non-MO Gallium Gadolinium Garnet (𝐺𝐺𝐺) substrate with 𝜀 𝐺𝐺𝐺 = 1.97 2 (see Fig. 1). The MO 𝐵𝐼𝐺 waveguide core is characterized by the complex dielectric tensor 𝜀 = 𝜀 ′ + 𝑖𝜀 ′′ . We assume that this material has MO activity due to a magnetization in the same plane as the ridge sides as noted in Fig. 1 (MT along the y direction). At optical frequencies, the MO activity in this configuration results in an optical anisotropy that can be described by the following dielectric tensor [START_REF] Vinovsky | Magneto-optical permittivity tensor in crystals[END_REF]:

𝜀 = ( 𝜀 𝑥𝑥 0 𝜀 𝑥𝑧 0 𝜀 𝑦𝑦 0 𝜀 𝑧𝑥 0 𝜀 𝑧𝑧 ) , (1)
Having considered a negligible absorption (𝜀 ′′ ≅ 0) for the BIG, the diagonal elements 𝜀 𝑥𝑥 , 𝜀 𝑦𝑦 , 𝜀 𝑧𝑧 are equal to the scalar isotropic dielectric constant 𝜀 𝐵𝐼𝐺 = 2.3 2 . The off-diagonal elements are 𝜀 𝑥𝑧 = -𝜀 𝑧𝑥 = 𝑖𝑔 where 𝑔 is the gyrotropy constant. 𝑔 is proportional to the magnitude of the magnetization. To better display the mechanisms undergoing in this structure, we first consider an artificial value of g as high as 0.1, assumed to be constant over the frequency range of interest. Such a high value of gyrotropy is chosen to highlight the device behavior and the performance trend. In the last part of the paper, we represent isolation ratio and FoM curves for different more realistic gyrotropy constants [START_REF] Vertruyen | Curie temperature exchange integrals, and magneto-optical properties in off-stoichiometric bismuth iron garnet epitaxial films[END_REF][START_REF] Doorman | Optical and magneto-optical tensor of bismuth-substituted yttrium iron-garnet film[END_REF].

The gold grating is adjacent to the 𝐵𝐼𝐺 waveguide. Each slit formed between the gold grating teeth acts as resonant MIM (Metal-Insulator-Metal) Fabry-Perot (FP) cavity which can enhance or suppress the transmission of light passing through it. In this analysis, we will consider grating slits ranging from 0.02 𝜇𝑚 to 0.06 𝜇𝑚 and grating period of 0.3 𝜇𝑚. The frequency dependent dielectric constant of gold is modeled by the Drude-Lorentz model fitting of ellipsometric data [START_REF] Palik | Handbook of optical constants of solids[END_REF] with plasma frequency 𝑤 𝑝 = 1.29 × 10 16 𝑟𝑎𝑑 𝑠 ⁄ .

In order to study the coupling of the gold grating to the MO waveguide and its role in achieving isolation, 3D simulations have been performed using the commercial software Lumerical FDTD Solutions. In the FDTD method, the isolator structure is described by a spatial discretization mesh. A spatial mesh size of 𝛥𝑥 = 𝛥𝑦 = 𝛥𝑧 = 10 𝑛𝑚 in the FDTD algorithm is found enough for convergence of the numerical results. We performed two kinds of simulations.

For the first type of simulations, we consider a single unit cell (one grating period including one resonant MIM structure) of the whole structure with periodic boundary conditions set in the propagation direction (see Fig. 2), 𝐵𝐼𝐺 and 𝑆𝑖𝑂 2 media are the z direction borders and account for the waveguide and surrounding medium respectively. In order to determine the waveguide modes present in the structure, as well as the possible resonant modes in the grating slits, we impose Bloch boundary conditions at the edge of the unitary cell along the propagation direction as represented in Fig. 2(a). This gives us the complete non-reciprocal band diagram of the infinite periodic structure.

Additionally, from such calculations, by properly inspecting the shape of the electromagnetic fields we can deduce the order of the resonant modes in the MIM structure as well as their nature. A second type of simulations is considered to investigate the full capabilities of this isolator, involving the whole structure with ten periodic gold gratings, placed adjacent to the waveguide and along the propagation direction. The waveguide is excited by a TE light source and PML (Perfectly Matched Layer) boundary conditions are applied at the boundaries of the simulation region to eliminate reflection of outgoing waves. This allows us to calculate the transmission spectra both in the forward and in the backward directions, and to deduce the isolation ratio and insertion losses. To obtain a high-performance isolator, we want (i) the insertion loss value to be as low as possible and (ii) the isolation ratio to be as high as possible [START_REF][END_REF]: the isolator figure of merit (FoM) is defined as the ratio of the former two.

Modes analysis

In the following three subsections, the propagative modes involved, their resonances and their impact on the behavior of the device are presented upon variation of the gold grating geometrical parameters, namely its period, its width and its length. Firstly, we evaluate the simplest geometry, i.e. by considering the isolator having closed slits (𝑙 𝑠 = 0): in this case the grating becomes a flat strip of gold supporting SPPs. Secondly, we examine the nature, the order and the spectral behavior of FP resonances occurring within the grating slits of finite size. Finally, the whole structure involving SPP-slit resonances coupling is analyzed.

SPP modes

The TE mode injected in the waveguide enables the excitation of SPPs at the interface of the waveguide-metal in the area of the grating. To study the effect of the gold on the waveguide mode, a simple structure made of a strip of gold placed adjacent to the waveguide is considered. The photonic band diagrams of a sole waveguide and a waveguide with strip of gold are presented in Fig. 3 (a) and (b) respectively.

Fig. 3: FDTD calculated band structures of (a) the sole dielectric waveguide and (b) the waveguide with a gold continuous strip (period 𝑝 𝐴𝑢 gold grating with 𝑙 𝑠 = 0). From top to bottom, the magnetic field profiles of (c) the fundamental 𝑇𝐸 𝑚𝑜𝑑𝑒 in case of (a), and of (d,e) the hybrid plasmonic-dielectric waveguide modes: 𝑀𝑜𝑑𝑒 𝐴 and 𝑀𝑜𝑑𝑒 𝐵.

Dash-dotted white lines correspond to 𝐵𝐼𝐺, 𝐺𝐺𝐺 and 𝑆𝑖𝑂 2 lightlines (from right to left). Fig. 3(a) shows the fundamental TE mode dispersion curve. It is worth pointing out that, due to the periodic boundary conditions used in the calculation, the dispersion curve is reflected at the edge of the First Brillouin Zone (FBZ). The upper picture in Fig. 3(c) shows the magnetic field profile of this fundamental TE mode.

The addition of a flat strip of gold next to the BIG waveguide allows the excitation of two hybrid plasmonic-dielectric modes (shown in Fig. 3(b)): the modes termed A (higher effective index) and B (lower effective index) display a more plasmonic and dielectric "character", respectively. Because of the lower index mismatch, the TE mode injected at the beginning of the waveguide will preferentially couple to the second hybrid waveguide plasmonic mode (𝑀𝑜𝑑𝑒 𝐵) at the vicinity of the gold strip. Moreover, the plasmonic nature of the both hybrid modes can enhance the TMOKE effect at the interface at any incident wavelength thanks to electromagnetic field concentration at the interface. It is worth pointing out the SPP is more concentrated with increasing frequency, when it approaches the plasma frequency (as a general property of SPPs).

Slit modes

Here we will consider the case of the grating having nontrivial slits (𝑙 𝑠 > 0). In the following FDTD calculations, dipole sources placed inside the 𝐵𝐼𝐺 waveguide can excite FP resonances within the slits of the gold grating. In order to determine these resonant modes in the considered frequency range, we will consider cavities (as presented in Fig. 2(b)) having a length 𝑙 𝑠 = 0.04 𝜇𝑚 and widths ranging between 0.1 and 0.9 𝜇𝑚 (z direction) and height equal to 0.45 𝜇𝑚 (y direction). A set of simulations is performed as previously explained and the resulting loci of the resonant frequencies are reported in Fig. 4. By inspecting Fig. 4, we can notice that, for a given frequency, the number of resonant peaks increases with the slit width (along the z direction). More generally, when the width and the height of the metal grating teeth are large enough the slit has the capability to host higher-order horizontal or vertical FP modes. Fig. 4 shows the evolution of the spectral position of these resonant modes with respect to the cavity width. These FP resonances are marked as 𝐻 𝑖 for horizontal resonance and 𝑉 𝑖 for modes including a vertical resonance respectively, where 𝑖 is related to the number of nodes in the magnetic field distributions [START_REF] Yoon | Analytic theory of the resonance properties of metallic nanoslit arrays[END_REF].

The spectral position of the horizontal cavity resonances in the z direction can be well fitted with the FP equation expressing the round-trip resonance when the cavity height is short enough to disable vertical high order modes and thus to avoid mode mixing [START_REF] Ding | Mapping surface plasmon polaritons and cavity modes in extraordinary optical transmission[END_REF]:

2𝑛𝜋 = 𝜑 1 + 𝜑 3 + 2𝑘 0 𝑛 𝑒𝑓𝑓,𝑠(𝑥,𝑦) 𝑤 𝐴𝑢, 𝑛 ∈ ℤ ( 2
)
where n is the order of the horizontal FP mode, 𝑤 𝐴𝑢 is the width of gold grating, 𝑛 𝑒𝑓𝑓,𝑠(𝑥,𝑦) is the effective index of the fundamental mode within the slits of the unitary cell calculated by considering a cross section parallel to the (x,y) plane of Fig. 2 and 𝜑 𝑟𝑖 is the reflection phase shift of the slit mode at each end of the cavity; 𝑖 = 1 and 𝑖 = 3 when the reflection of light from the 𝐵𝐼𝐺 waveguide and from the 𝑆𝑖𝑂 2 superstrate are considered, respectively. In the same manner, the spectral position of vertical FP resonances can be modelled accurately when the cavity width is short enough to host only vertical high order resonances. In this case the equation is as follows:

2𝑚𝜋 = 𝜑 2 + 𝜑 3 + 2𝑘 0 𝑛 𝑒𝑓𝑓,𝑠(𝑥,𝑧) ℎ 𝐴𝑢, 𝑚 ∈ ℤ (3)
where m is the order of the vertical FP mode, ℎ 𝐴𝑢 is the height of gold grating, 𝑛 𝑒𝑓𝑓,𝑠(𝑥,𝑧) is the effective index of the fundamental mode within the slits of the unitary cell calculated by considering a cross section parallel to the (x,z) plane of Fig. 2 and 𝜑 𝑟𝑖 is the reflection phase shift of the slit mode at each end of the cavity; 𝑖=2 when the reflection from the 𝐺𝐺𝐺 substrate is considered. It is worth pointing out that the reflection phase shift can be retrieved from the argument of the complex coefficient 𝑟 𝑖 , which at each boundary can be roughly approximated as the normal incidence reflection coefficient.

𝑟 𝑖 = |𝑟 𝑖 | exp(𝑖𝜑 𝑟𝑖 ) = 𝑛 𝑒𝑓𝑓,𝑠 -√𝜀 𝑖 𝑛 𝑒𝑓𝑓,𝑠 + √𝜀 𝑖 , (4)
For the narrowest cavity width (𝑤 𝐴𝑢 = 0.1 𝜇𝑚) and in the frequency range of interest, we can notice a single resonant mode which is a fundamental order vertical FP mode 𝑉 0 . Increasing the cavity width to 𝑤 𝐴𝑢 = 0.3 𝜇𝑚 introduces a fundamental order horizontal FP 𝐻 0 and a first order vertical FP 𝑉 1 resonance, at 155 THz and 195 THz respectively. With further increase of the cavity width from 𝑤 𝐴𝑢 = 0.3 𝜇𝑚 up to the highest width 𝑤 𝐴𝑢 = 0.9 𝜇𝑚, six resonant modes finally appear. After inspection of the magnetic field profiles inside the slits, these modes have been identified as the 𝐻 0 , 𝑉 0 , 𝐻 1 , 𝑉 1 , 𝑉 2 and 𝐻 2 FP resonances. Nevertheless, in the range of calculation of Fig. 4, the 1 st and 2 nd order "vertical" FP modes 𝑉 1 and 𝑉 2 undergo a distortion of the profile as shown in the subplots (c) and (d). Whatever the mode shape, the response of the FP resonant cavity can be controlled by tuning the geometrical parameters of the grating: this can be fruitfully exploited to control the frequency at which isolation will occur for the whole structure.

SPP-slit modes coupling

When the side gold grating of the full structure (Fig. 2) has nontrivial slits (𝑙 𝑠 > 0), the hybrid plasmonic-dielectric mode (𝑀𝑜𝑑𝑒 𝐴 or 𝑀𝑜𝑑𝑒 𝐵 of Fig. 3) couples to the slits FP resonances. Similarly to the non-guided version [START_REF] Halagacka | Magnetoplasmonic nanograting enables optical nonreciprocity sign reversal[END_REF], such coupling will lead to the enhancement of TMOKE effect and thus of propagation non-reciprocity.

To shed light on the coupling mechanisms between the waveguide and the plasmonic grating, the photonic band diagram of the complete infinite periodic structure (grating and waveguide as depicted in Band diagrams shown in Fig. 5(a) and (c) are representative of the logarithm of the excited modes intensity. In the corresponding handmade tracing reported in Fig. 5(b) and (d), the red and the blue dispersion curves correspond to the magneto-plasmonic waveguide modes A and B, respectively, whereas the "flat" curves represent different orders of FP resonance modes. The gap at 230 THz at the FBZ border corresponds to the first order Bragg grating reflection, which depends here mainly on the grating period (0.3 𝜇𝑚 in the both cases). Since reflections must be avoided in isolators, we will focus our analysis on shorter frequency ranges. Besides, below 150 THz we know from Fig. 2(a,b) that the guided mode is mainly plasmonic (SPP-like 𝑀𝑜𝑑𝑒 𝐴). Such mode can efficiently interact with FP slit modes: this is manifested by the dispersion curve bending of 𝐻 0 , 𝑉 0 , 𝐻 1 and 𝑉 1 which attest of interaction with weakly confined or leaky waveguide modes. Nevertheless, the frequency of most interest here is between 150 THz and 220 THz, i.e. the NIR range. In this range, one clear anticrossing is identified in the case 𝑤 𝐴𝑢 = 0.1 𝜇𝑚 (Fig. 5(a)) between the waveguide 𝑀𝑜𝑑𝑒 𝐵 (mainly) and the slit mode 𝑉 1 (blue and green curves of the schematic view). Similarly, a gap is opened near 180 THz in the 𝑤 𝐴𝑢 = 0.9 𝜇𝑚 case (Fig. 5(b)), resulting from the interaction between the waveguide 𝑀𝑜𝑑𝑒 𝐴 (mainly) and the slit FP resonance 𝐻 2 (burgundy and red curves of the schematic view).

As shown in the next paragraph, the intersection of these bands can lead to resonant dips in the transmission spectra, without inducing reflection, which are in turn the points at which enhanced isolation occurs.

Isolator analysis

The full guided structure with a ten periods gold grating located on the side of the 𝐵𝐼𝐺 waveguide has been simulated for two cases: 𝑤 𝐴𝑢 = 0.1 𝜇𝑚 and 𝑤 𝐴𝑢 = 0.9 𝜇𝑚. The grating is sized with the same parameters as stated before (period of 0.3 𝜇𝑚) and has thus a total length of 3 𝜇𝑚 . Forward and backward transmissions 𝑇 𝑓 and 𝑇 𝑏 are calculated as well as the isolation ratio, defined as:

𝐼𝑅 (𝑑𝐵) = 10 log 10 ( 𝑇 𝑓 𝑇 𝑏 ⁄ ), (5)
Fig. 6(a) represents the isolation spectra obtained in the both cases. Each curve contains one or several S-shaped peaks which demonstrate the enhancement of the non-reciprocal behavior at given frequencies. Interestingly the nonreciprocity sign can be positive or negative, despite the same magnetization orientation. Besides the number of Sshaped peaks increases with the width of the cavity. These S-shaped peaks are similar to the ones described in the non-guided structure reported in [START_REF] Halagacka | Magnetoplasmonic nanograting enables optical nonreciprocity sign reversal[END_REF] and can be attributed to the magneto-optical interaction enhancement by the SPP-slit mode coupling. Their increased number in the wider cavity is due to the existence of higher order FP resonances hosted in the grating slits, as explained in the paragraph 3.2. When 𝑤 𝐴𝑢 = 0.9 𝜇𝑚, the isolation ratio curve (in black) shows three prominent resonant S-shaped features labeled 1,2 and 3 at 125 THz, 185 THz and 245 THz, respectively. By looking at the magnetic field of the full structure shown in Fig. 7, we can clearly identify that these points correspond to the 1 st horizontal order (a), the 2 nd horizontal order (b) and the 3 rd horizontal order (c) Fabry-Perot modes. A closer look at the band diagrams of Fig. 5 allows us better understanding the amplitude difference at the different isolation points.

In the case of 𝑤 𝐴𝑢 = 0.1 𝜇𝑚, the interaction of the 1 st order vertical FP mode occurs mainly with the hybrid waveguide 𝑀𝑜𝑑𝑒 𝐵 (marked by 𝑀𝑜𝑑𝑒 𝐵 and 𝑉 1 in Fig. 5(a)) and is manifested by the anticrossing near 220 THz. The coupling between the hybrid waveguide mode and the slit FP resonance results in transmission non-reciprocity enhancement as seen on the red curve of Fig. 6(a). The magnetic field intensity of the FP resonance occurring in the grating slits at this high isolation point is shown in Fig. 8. We have checked that the Bragg reflection (which occurs at 230 THz) is negligible at this frequency, and thus optical isolation is fully performed. Since the SPP is more concentrated at higher frequency, this is a particularly favorable configuration combining enhanced TMOKE, high confinement of (short slit width) FP mode and sub-wavelength grating period. In the case of 𝑤 𝐴𝑢 = 0.9 𝜇𝑚 (in black), three peaks appear on the isolation curve. Another important parameter which quantitatively describes the efficiency of the isolator is the FoM defined as the ratio of isolation to the insertion losses:

𝐹𝑜𝑀 = 𝐼𝑅 (𝑑𝐵) -10 log 10 (𝑇 𝑓 ) , (6)
FoM calculations are shown in Fig. 6(b) for both cavity widths. It is worth pointing out that the FoM decreases with increasing cavity widths. This can be attributed to losses induced in larger cavities FP modes, which do not contribute directly to non-reciprocal effect. Besides the highest FoM obtained for 𝑤 𝐴𝑢 = 0.1 𝜇𝑚 at 218.5 THz, when 𝑔 = 0.1 is considered, corresponds to 22 𝑑𝐵 isolation ratio, 17.42 𝑑𝐵 insertion loss and negligible back reflection. In addition to the higher compactness of the structure, the grating isolator with the shortest cavity width 𝑤 𝐴𝑢 = 0.1 𝜇𝑚 presents promising performance.

Furthermore, we intend to optimize both the isolation ratio and the FoM by varying the grating duty cycle. The period 𝑝 𝐴𝑢 remains fixed at 0.3 𝜇𝑚, and the slit length 𝑙 𝑠 is varied. The positions of the isolation peaks undergo a blueshift as the slit lengths increase from 0.02 𝜇𝑚 to 0.06 𝜇𝑚 as shown in Fig. 9(a), for a gold grating with ten periods and wAu= 0.1 𝜇𝑚. Since the isolation ratio enhancement results from the hybrid plasmonic modes and the FP slit mode coupling, the amplitude balance between the coupled modes is reproduced on the isolation ratio peaks when the FP mode frequency shifts. Thus, the positive and the negative non-reciprocity may have different absolute values. Besides in this range, the FoMs reach similar acceptable levels (>1). It appears here that the variation of the slit length induces a slight isolation frequency shift. Besides, we have shown that isolation can be broadly tuned by the slit width (Fig. 4 and6). Thus, the frequency at which isolation occurs can be coarsely tuned by playing on slit width and finely tuned via slight adjustment of the grating slit length. In this way, grating geometry allows covering a wide and continuous range of isolation frequency, with similar FoM as long as the FP slit modes appear in the frequency range of waveguide hybrid modes and SPP mode coupling.

Isolator performance evaluation

In the previous parts we explained the intertwined mechanisms involved in providing the isolation effect in the proposed magnetoplasmonic structure. We have also assessed the influence of the geometrical parameters, by using an artificially high and dispersionless gyrotropy value. Although this can depict the performance trend for high gyrotropy, here we consider the effect of more usual gyrotropy values [START_REF] Srinivasan | High-Gyrotropy Seedlayer-Free Ce:TbIG for Monolithic Laser Matched SOI Optical Isolators[END_REF][START_REF] Bsawmaii | Longitudinal magneto-optical effect enhancement with high transmission through a 1D all-dielectric resonant guided mode grating[END_REF] on the isolation ratio and FoM curves.

For the sake of simplicity, we consider a constant value on the frequency range of interest. Fig. 10 shows the isolation ratio as well as the FoM magnitudes for different gyrotropy levels. Lowering the gyrotropy shrinks the bandwidth of isolation, from 15.26 𝑛𝑚 for 𝑔 = 0.1 to 2.23 𝑛𝑚 for 𝑔 = 0.01 and 1.95 𝑛𝑚 for 𝑔 = 0.005. The bandwidth shrinking is due to the decrease in the shift of the FP resonance position in the opposite senses. Although the isolation mechanism presented in this manuscript does not depend on the gyrotropy value, lower value of gyrotropy reduces the device performance by decreasing both isolation and bandwidth by effectively decreasing the variation between the propagation constants in the forward and backward propagation directions due to the TMOKE effect. Nevertheless, with no specific optimization, significant isolation ratio higher than 5 𝑑𝐵 remains reachable for 𝑔 = 0.005 with a grating as short as 3 𝜇𝑚. A larger number of grating elements combined with overall optimization of the geometrical parameters such as grating width and slit length could re-enhance the isolation and increase its ratio.

Conclusion

In conclusion, we have proposed a new isolator design for TE modes based on magnetoplasmonic coupling between a magnetooptical waveguide and a side gold subwavelength grating. This very compact structure, with a total length of 3 𝜇𝑚, can reach typically from 5 𝑑𝐵 to 22 𝑑𝐵 isolation ratio with a gyrotropy from 0.005 to 0.1, and negligible back reflection at telecom wavelength. Fabry-Perot resonances in the grating slits coupled to the SPP component of the waveguide mode lead to this strong enhancement of the isolation ratio. A change in the grating geometrical parameters enables the selection of the wavelength of interest by concurrently granting significant isolation on a wide wavelength range. Besides, the isolation ratio can be positive or negative with the same magnetization orientation, enabling integration of isolators in different propagation directions on the same circuit with a single external magnetic field. In summary, this very compact and waveguiding structure may pave the way for achieving novel compact optical isolators and can simultaneously provide a versatile solution for the integration of nonreciprocal components in photonic circuits.

Fig. 1 :

 1 Fig. 1: Scheme of the proposed magneto-plasmonic isolator. A periodic gold grating is positioned on one side of a 𝐵𝐼𝐺 waveguide integrated on a non-magnetic 𝐺𝐺𝐺 substrate. MT is the magnetization of the BIG layer, oriented along the y direction.

Fig. 2 :

 2 Fig. 2: Sketch of the device unitary cell. (a) Axonometric view of the unitary cell. (b, c, d) 2D cut views. In (b) the upper and lower horizonal dashed thin lines represent the traces of the planes where the 2D cuts depicted in (c) and in (d) are taken, respectively. The two vertical thin dotted lines represent the position where the Bloch boundary conditions are imposed in FDTD calculation to delimit the unitary cell. In (c) and (d) the horizonal dashed thin line represents the trace of the plane where the 2D cut depicted in (b) is taken.

Fig. 4 :

 4 Fig.4: (a) Dispersion curves of the different horizontal and vertical cavity modes FP resonances (marked as 𝐻 𝑖 or 𝑉 𝑖 , respectively) within the grating slits as a function of the cavity width (z direction). (b)-(g) Intensity of the magnetic field along the (y,z) plane, when x corresponds to the center of the slits, for the different points marked with black circles in (a).

Fig 2 )

 2 has been calculated. Fig. 5(a-d) represent the band diagrams when the grating height and the grating period are fixed at ℎ 𝐴𝑢 = 0.45 𝜇𝑚 and 𝑝 𝐴𝑢 = 0.3 𝜇𝑚, the grating length equals 𝑙 𝑠 = 0.04 𝜇𝑚, and the widths are (a,b) 𝑤 𝐴𝑢 = 0.1 𝜇𝑚 and (c,d) 𝑤 𝐴𝑢 = 0.9 𝜇𝑚 respectively.

Fig. 5 :

 5 Fig. 5: Simulated photonic band diagrams for (a) 𝑤 𝐴𝑢 = 0.1 𝜇𝑚 and (c) 𝑤 𝐴𝑢 = 0.9 𝜇𝑚. To improve their readability, a handmade tracing of the corresponding dispersion curves are reported in (b) and (d) respectively.

Fig. 6 :

 6 Fig. 6: (a) Isolation ratio and (b) FoM for 𝑤 𝐴𝑢 = 0.1 𝜇𝑚 (in red) and 𝑤 𝐴𝑢 = 0.9 𝜇𝑚 (in black), when ℎ 𝐴𝑢 = 0.45 𝜇𝑚.

Fig. 7 :

 7 Fig. 7: First row: (a-c) Intensity of the magnetic field along the (x,z) plane (top view of the device), when 𝑦 = 0 𝜇𝑚. Second row: (d-e) magnetic field intensity maps along the (y,z) plane, when 𝑥 = 0 𝜇𝑚. Each column corresponds respectively to the three isolation frequencies marked on the black curve of Fig. 6 as 1, 2 and 3: at (a) 125 THz, (b) 185 THz and (c) 245 THz which involve modes 𝐻 1 , 𝐻 2 and 𝐻 3 in Fig. 5(b).

  Points 1, 2 and 3 are related to 𝐻 1 , 𝐻 2 and 𝐻 3 resonances respectively (represented in Fig.8(a,b,c)). By analyzing the band diagram (Fig.5(c)), 𝐻 1 and 𝐻 3 resonances appear weakly coupled to the waveguide mode without causing any band gap, especially because the waveguide mode is very deconfined and leaky below 150 THz and undergoes reflections above 230 THz. The interesting point is at 185 THz where a 2 nd order horizontal FP (𝐻 2 ) intersects the two hybrid Modes A and B dispersion curves near their anticrossing gap: at this frequency, the coupling of the dielectric waveguide mode and the SPP mode is already maximum (Fig.3(b)), and thus favors their interaction, and especially the interaction of Mode A, with the FP slit mode. This could explain the increase in isolation ratio near 185 THz. However, the maximum isolation ratio is lower than in the case of 𝑤 𝐴𝑢 = 0.1 𝜇𝑚, probably because of less confined SPP at shortest wavelength, and less confined FP mode in the wider cavity.

Fig. 8 :

 8 Fig. 8: Intensity of the magnetic field (a) along the (x,z) plane (top view of the device), at 𝑦 = 0 𝜇𝑚 and (b) along the (x,y) plane (side view of the device), at 𝑧 = 0.55 𝜇𝑚, calculated at 220 THz. ℎ 𝐴𝑢 = 0.45 𝜇𝑚 and 𝑤 𝐴𝑢 = 0.1 𝜇𝑚.

Fig. 9 :

 9 Fig. 9: (a) Isolation ratio in dB (color scale) (b) FoM for 𝑤 𝐴𝑢 = 0.1 𝜇𝑚 and 𝑝 𝐴𝑢 = 0.3 𝜇𝑚 versus slit length for a device length of 3 𝜇𝑚.

Fig. 10 :

 10 Fig. 10: (a) Isolation ratio in dB (color scale) (b) FoM for 𝑤 𝐴𝑢 = 0.1 𝜇𝑚 and 𝑝 𝐴𝑢 = 0.3 𝜇𝑚 for gyrotropy values 𝑔 = 0.1, 0.01, 𝑎𝑛𝑑 0.005.

Acknowledgements

No potential conflict of interest was reported by the authors.