Franck Polewczyk

Jean-Marc Leyssale

Paul Lafourcade
email: paul.lafourcade@cea.fr

Temperature-dependent elasticity of single crystalline graphite

Despite its relevance in many high temperature processes, the elastic behavior of single crystalline graphite is so far entirely undocumented away from room temperature conditions. In this work we present a molecular dynamics investigation of graphite's second order elastic tensor dependence on temperature from 300 to 4000 K, using a series of popular interatomic potentials for carbon. Data and analytic expressions are presented for the elastic tensor under the two limiting situations known as isothermal, or quasi-static, and adiabatic deformation conditions, prevailing in the limits of slow and large deformation rates, respectively. Independently from the potential, we identify a strong non-linearity of elastic constants with respect to temperature. We show that despite conserving an important elastic anisotropy whatever the temperature, the latter is being reduced by a factor of ∼ 5 when increasing temperature up to 4000 K. Also, we show that elastic anisotropy is about 20 times larger under isothermal conditions than under adiabatic conditions. Finally, we investigate the dynamics of the adiabatic to isothermal stress relaxations occurring right after ultra fast deformations such as those encountered under shock loading, showing that the relevant time and length scales are of the order of a few ps and nm, respectively, for deformations around a percent.

Introduction

When considering structural elements for devices experiencing ultra-high temperatures (UHT), like for instance thermal protection systems for atmospheric re-entry [START_REF] Couzi | Improvements in Ablation Predictions for Reentry-Vehicles Nosetip[END_REF], involving temperatures in the 3000-4000 K range, very few materials are available [START_REF] Shabalin | Ultra-high temperature materials II[END_REF]. Carbon based materials, and especially carbon/carbon (C/C) composites, combining high melting temperature, stiffness, strength and thermal conductivity, as well as low weight, are thus generally selected [START_REF] Savage | Carbon/Carbon composites[END_REF][START_REF] Scarponi | Carbon-carbon composites in aerospace engineering[END_REF]. However, despite numerous investigations, our knowledge of the properties of C/C composites and of their constituents, the well-known carbon fibers and pyrolytic carbon (pyC) matrices, remains limited to significantly lower temperatures. [START_REF] Trinquecoste | High temperature thermal and mechanical properties of high tensile carbon single filaments[END_REF][START_REF] Sauder | Thermomechanical properties of carbon fibres at high temperatures (up to 2000 °C)[END_REF][START_REF] Sauder | The tensile behavior of carbon fibers at high temperatures up to 2400[END_REF][START_REF] Sauder | The tensile properties of carbon matrices at temperatures up to[END_REF][START_REF] Pradere | Transverse and longitudinal coefficient of thermal expansion of carbon fibers at high temperatures (300-2500k[END_REF]. Therefore, the properties of interest -coefficients of thermal expansion, Young's moduli and Poisson coefficients along and normal to the fiber axes, etc. -have to be extrapolated to the temperatures of interest, which can be somehow heuristic.

Due to these extreme conditions that can be challenging, experiments are sometimes not suitable to identify and understand complex behaviors involved at the atomic and mesoscopic scales. To overcome this, virtual material approaches are becoming more and more popular.

These approaches have indeed proven their value in understanding the thermo-mechanical and chemical behavior of high explosives [START_REF] Kroonblawd | High explosive ignition through chemically activated nanoscale shear bands[END_REF][START_REF] Kroonblawd | Anisotropic strength behavior of single-crystal tatb[END_REF][START_REF] Lafourcade | Mesoscopic constitutive law with nonlinear elasticity and phase transformation for the twinning-buckling of tatb under dynamic loading[END_REF], in investigating metals response under very high strain rates [START_REF] Zepeda-Ruiz | Probing the limits of metal plasticity with molecular dynamics simulations[END_REF][START_REF] Zepeda-Ruiz | Atomistic insights into metal hardening[END_REF][START_REF] Van Der Giessen | Roadmap on multiscale materials modeling[END_REF], as well as in studying the damage and mechanical behavior of ceramics and carbon composites [START_REF] Mazars | Damage investigation and modeling of 3d woven ceramic matrix composites from x-ray tomography in-situ tensile tests[END_REF][START_REF] Gillard | Modeling of the non-linear mechanical and thermomechanical behavior of 3d carbon/carbon composites based on internal interfaces[END_REF].

When bridging the gap between time and length scales, which is the essence of multiscale modeling, complex material architectures -as those in C/C composites -can be simulated under various thermodynamic conditions, provided one has a precise enough constitutive law for the individual elements of the microstructure. However, determining such constitutive laws under UHT conditions can be even more challenging than investigating the whole composite material, from the experimental point of view.

Atomistic simulation techniques, like molecular dynamics (MD) simulations, are, in theory, well-suited for such investigations. In the case of C/C composites, large-scale (2 × 10 5 atoms) atomistic models of the pyC matrices [START_REF] Farbos | Nanoscale structure and texture of highly anisotropic pyrocarbons revisited with transmission electron microscopy, image processing, neutron diffraction and atomistic modelling[END_REF] have been proposed in recent years. However, their corresponding room temperature elastic constants [START_REF] Farbos | Nanoscale elasticity of highly anisotropic pyrocarbons[END_REF] failed to reproduce the existing experimental data [START_REF] Sauder | The tensile properties of carbon matrices at temperatures up to[END_REF]. In particular, the longitudinal tensile moduli were significantly larger, by factors of 2 to 6, than the measurements, indicating that the models are either not accurate nor sufficiently large to properly capture the elastic behavior of the pyC. We note that, while currently no such atomistic model exist to describe the bulk properties of C fibers, some models of carbon fibers surfaces have been recently proposed [START_REF] Vuković | Practical atomistic models of carbon fiber surfaces with tuneable topology and topography[END_REF].

Taking some distance with the accurate representation of actual C/C composites, we consider in this work the case of crystalline hexagonal graphite, which, conversely to the case of pyC matrices and carbon fibers, has a perfectly determined structure. Also, experimental data only exists for the full elastic tensor of single crystalline graphite at room temperature [START_REF] Blakslee | Elastic constants of compressionannealed pyrolytic graphite[END_REF][START_REF] Bosak | Elasticity of single-crystalline graphite: Inelastic x-ray scattering study[END_REF], which is not the case for fibers and matrices for which only longitudinal and transverse Young's moduli are generally reported. On the other hand, to the best of our knowledge, no experimental data is available concerning the elastic constants of graphite at elevated temperatures.

In this work we use MD simulations to determine the temperature dependence of graphite's unit cell lengths and elastic constants, which are fundamental parameters in deriving constitutive laws for up-scaled thermo-mechanical models of graphite-based complex architectures, inaccessible to atomic-scale modeling. Considering applications under somehow extreme conditions such as shock loading, high strain-rate deformation or thermal aggression, these properties are computed up to ultra high temperature, i.e. T = 4000 K, at ambient pressure, i.e. P = 1 atm. Also, two types of elastic constants, corresponding to drastically different conditions, are determined. On one hand, we compute the usual isothermal (or quasi-static) elastic constants, in which strain rate is low with respect to all elastic relaxation modes of the material. On the other hand, the so-called adiabatic elastic constants are computed, under the application of instantaneous deformations without any subsequent atomic positions relaxation. This method has already been used to compute the adiabatic elastic constants of the triclinic compound TATB [START_REF] Lafourcade | Dislocation core structure at finite temperature inferred by molecular dynamics simulations for 1,3,5-triamino-2,4,6-trinitrobenzene single crystal[END_REF]. As explained by Wallace [START_REF] Wallace | Thermodynamics of Crystals[END_REF] and Sutton [START_REF] Sutton | Physics of Elasticity and Crystal Defects, Oxford Series on Materials Modelling[END_REF], any given thermoelastic process should be associated with either adiabatic or isothermal elastic constants. However, whether it should be one or the other is directly related to the rate of elastic deformation. For example, adiabatic elastic constants are measured by ultrasonic pulse experiments while their isothermal counterpart are usually measured in direct quasi-static mechanical tests at constant temperature. The difference between adiabatic and isothermal elastic constants is in general temperature dependent, which will be shown in the present work. In the end, analytical formulations of both elastic tensors will be provided as a function of temperature, ready to be assigned to a continuum framework. However, if the later involves a physical process that requires both adiabatic and isothermal elastic moduli over various thermodynamic conditions, it should always be associated with corresponding thermal stress/strain tensors as well as specific heat.

While extended benchmarks of carbon potentials have been published in recent years [START_REF] De Tomas | Graphitization of amorphous carbons: A comparative study of interatomic potentials[END_REF][START_REF] De Tomas | Transferability in interatomic potentials for carbon[END_REF], the latter focused on the structure of quenched, then annealed, disordered carbons. So far, the abilities of such potentials to describe the lattice parameters and elastic constants of graphite, and their evolution upon increasing temperature remains mostly unknown. Therefore, and as no experimental data exists to validate the simulations, we consider five popular empirical potentials: REBO [START_REF] Brenner | A second-generation reactive empirical bond order (rebo) potential energy expression for hydrocarbons[END_REF], AIREBO [START_REF] Stuart | A reactive potential for hydrocarbons with intermolecular interactions[END_REF], AIREBO-M [START_REF] O'connor | Airebo-m: A reactive model for hydrocarbons at extreme pressures[END_REF], CEDIP [START_REF] Marks | Generalizing the environment-dependent interaction potential for carbon[END_REF] and LCBOPII [START_REF] Los | Improved long-range reactive bond-order potential for carbon. i. construction[END_REF]. Comparison of the results obtained with these different potentials allows us to determine reliable trends on the high temperature elastic behavior of graphite.

Methods

Interatomic potentials

Empirical potentials are frequently used to investigate elastic and fracture properties of carbon-based systems [START_REF] Grantab | Anomalous strength characteristics of tilt grain boundaries in graphene[END_REF][START_REF] Zhang | Fracture toughness of graphene[END_REF][START_REF] Gamboa | Mechanism of strength reduction along the graphenization pathway[END_REF]. Here we consider five potentials amongst the most widely used: REBO [START_REF] Brenner | A second-generation reactive empirical bond order (rebo) potential energy expression for hydrocarbons[END_REF], AIREBO [START_REF] Stuart | A reactive potential for hydrocarbons with intermolecular interactions[END_REF], AIREBO-M [START_REF] O'connor | Airebo-m: A reactive model for hydrocarbons at extreme pressures[END_REF], CEDIP [START_REF] Marks | Generalizing the environment-dependent interaction potential for carbon[END_REF] and LCBOPII [START_REF] Los | Improved long-range reactive bond-order potential for carbon. i. construction[END_REF]. The reactive empirical bond order (REBO) potential of Brenner et al. [START_REF] Brenner | A second-generation reactive empirical bond order (rebo) potential energy expression for hydrocarbons[END_REF] has been extensively used in the literature, mostly to investigate elasticity and fracture of nanocarbons, including graphene and carbon nanotubes [START_REF] Zhang | Fracture toughness of graphene[END_REF][START_REF] Gamboa | On the prediction of graphene's elastic properties with reactive empirical bond order potentials[END_REF][START_REF] Lu | Atomistic simulation and continuum modeling of graphene nanoribbons under uniaxial tension[END_REF]. The two versions of the adaptive intermolecular reactive bond order potential (AIREBO and AIREBO-M) from Stuart et al. [START_REF] Stuart | A reactive potential for hydrocarbons with intermolecular interactions[END_REF] are directly derived from REBO. They complement the latter with a bond-order dependent integration of van der Waals interactions and a torsional potential around single C-C bonds to better describe hydrocarbons in the liquid state and aromatic stacking distances. In the AIREBO-M potential, the Lennard-Jones term used to describe van der Waals interactions in AIREBOknown to significantly overestimate repulsive forces under large compressive strains -is replaced with a more realistic Morse potential [START_REF] O'connor | Airebo-m: A reactive model for hydrocarbons at extreme pressures[END_REF]. Another potential, namely SED-REBO [START_REF] Perriot | Screened environment-dependent reactive empirical bond-order potential for atomistic simulations of carbon materials[END_REF] is available in the literature and based on the second generation of REBO. Since this potential is equivalent to REBO in terms of elasticity, it was not used in the present work. The last two potentials considered in the present work are the environment dependent interatomic potential for carbon (CEDIP) from Marks et al. [START_REF] Marks | Generalizing the environment-dependent interaction potential for carbon[END_REF] and the improved long-range carbon bond-order potential (LCBOPII) from Los et al. [START_REF] Los | Improved long-range reactive bond-order potential for carbon. i. construction[END_REF]. Conversely to the previously mentioned potentials in which bonding interactions are cut off at very short distances, via a switching function operating in the 1.7-2 Å range, medium-range interactions are included in the both CEDIP and LCBOPII potentials, which is a clear advantage for capturing fracture properties.

Another major consequence of medium-range interactions is that these potentials allow for a quantitative prediction of the energy barrier for the graphite to diamond transition while potentials of the REBO family do not [START_REF] Los | Formation of multiwall fullerenes from nanodiamonds studied by atomistic simulations[END_REF]. However, LCBOPII differs from CEDIP since it contains explicit long-range (dispersion) interactions while the former does not. Although these two potentials have been considerably less used than REBO or AIREBO potentials, a few reports suggest that these potentials are valuable candidates. Elastic and fracture properties of low temperature and annealed amorphous carbon have been investigated using CEDIP [START_REF] De Tomas | Carbide-derived carbons for dense and tunable 3d graphene networks[END_REF], while LCBOPII was used to study the pressure-induced graphite to diamond phase transformation under shock conditions and homothetic strain [START_REF] Pineau | Molecular dynamics simulations of shock compressed graphite[END_REF][START_REF] Pineau | Molecular dynamics simulations of shock compressed heterogeneous materials. ii. the graphite/diamond transition case for astrophysics applications[END_REF][START_REF] Lafourcade | Elastic instability in graphite single crystal under dynamic triaxial compression: Effect of strain-rate on the resulting microstructure[END_REF].

Molecular dynamics simulations

Two different classical molecular dynamics (MD) codes are used in this work. The LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [START_REF] Plimpton | Fast parallel algotithms for short-range molecular dynamics[END_REF][START_REF] Thompson | LAMMPS -a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales[END_REF] is used for all simulations using REBO, AIREBO, AIREBO-M and CEDIP potentials while the STAMP [START_REF] Soulard | Molecular Dynamics Study of the Micro-spallation[END_REF] (Simulations Temporelles Atomistiques Massivement Parallélisées) code is used to perform simulations using LCBOPII potential as it is not currently implemented in LAMMPS. For both codes, MD trajectories are integrated using a velocity-Verlet integrator [START_REF] Verlet | Computer" experiments" on classical fluids. ii. equilibrium correlation functions[END_REF] and a 0.5 fs timestep, ensuring a correct energy conservation under equilibrium microcanonical (NVE) runs. Besides energy conservation tests, simulations are run in both canonical (NVT) and isothermal-isobaric (NPT) ensembles, using a Nosé-Hoover style thermostat [START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Hoover | Canonical dynamics: Equilibrium phase-space distributions[END_REF] and barostat [START_REF] Martyna | Constant pressure molecular dynamics algorithms[END_REF] along with 0.1 ps and 1.0 ps damping constants, respectively. For the determination of graphite single crystal lattice parameters and elastic constants evolution with temperature, a 3D-periodic orthorhombic simulation cell containing 3456 carbon atoms is considered, consisting in a 6 × 6 × 6 replication of the 16 atoms graphite orthorhombic unit cell. In Section 3.6, a 107520 atoms 3D-periodic simulation cell consisting in a 20 × 24 × 14 replication of the same unit cell is considered to investigate the dynamical aspects of stress relaxation with sufficient statistics. In all the simulations, except the ones for computing the elastic constants that require a change in the simulation box shape, i.e. when applying shear strain, the cell angles are constrained to maintain a constant value of 90 %, accordingly to the graphite orthorhombic unit cell. For the simulations in the isothermal-isobaric ensemble, diagonal components of the stress tensor are controlled independently so that true hydrostatic conditions can be achieved, accordingly to the very large elastic anisotropy of graphite single crystal.

Lattice equilibration at finite temperature

For each potential, the system is equilibrated in the isothermal-isobaric (NPT) ensemble during 1 ns at a pressure of 1 atm and at 300 K and for temperatures ranging from 500 to 4000 K (i.e. close to graphite's melting point), by steps of 500 K. Lattice parameters are systematically averaged over the last 250 ps of the NPT trajectory, during which the stress tensor is checked to be fully hydrostatic, without any substantial evolution of simulation cell lengths. As both REBO and CEDIP potentials do not include van der Waals interactions, only the in-plane lattice parameters are allowed to evolve in the NPT simulations for these two potentials, the out-of-plane parameter were arbitrarily held fixed to the value obtained using the LCBOPII model at the considered temperature.

Elastic constants calculation

Graphite single crystal elasticity is known to possess a transverse isotropy symmetry, which is an elastic symmetry derived from the orthotropic elastic symmetry, for which the second-order elastic tensor possesses 9 independent elastic components:

C =                               C 11 C 12 C 13 0 0 0 C 12 C 22 C 23 0 0 0 C 13 C 23 C 33 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 55 0 0 0 0 0 0 C 66                               (1)
where the Voigt contraction convention has been used to represent the elastic tensor C as a second-order tensor. As an additional criterion to compare the results obtained with the different potentials, we first make the assumption that the graphite single crystal has an orthotropic symmetry, with 9 independent elastic constants. This hypothesis is always true because a material that possesses a transverse isotropic elasticity has by definition an orthotropic elasticity. However, for a material to have a transverse isotropic elasticity, it is

required that, C 11 = C 22 , C 13 = C 23 , C 44 = C 55 and C 66 = (C 11 -C 12 )/2
, with all other components equal to zero. For graphite single crystal this implies that the elastic response to a deformation along the armchair or zigzag direction is strictly equivalent, which might not be true for certain potentials [START_REF] Gamboa | On the prediction of graphene's elastic properties with reactive empirical bond order potentials[END_REF]. When a material is transversely isotropic, its second-order elastic tensor reads:

C =                               C 11 C 12 C 13 0 0 0 C 12 C 11 C 13 0 0 0 C 13 C 13 C 33 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 44 0 0 0 0 0 0 (C 11 -C 12 )/2

Isothermal elastic constants

Isothermal elastic constants are computed as follows. Starting from the equilibrium lattice parameters, a series of uniaxial engineering strains are uniformly applied to the graphite supercell. The considered strain range correspond to ε ∈ [-0.3, +1.0] %, with a 0.1 % increment, for deformations along the principal axes, where negative and positive values correspond to compressive and tensile strains, respectively. Note that the lower range of compressive strain was chosen to avoid well-known compression-induced buckling elastic instabilities in anisotropic materials [START_REF] Lafourcade | Irreversible deformation mechanisms for 1,3,5-triamino-2,4,6trinitrobenzene single crystal through molecular dynamics simulations[END_REF][START_REF] Lafourcade | Mesoscopic constitutive law with nonlinear elasticity and phase transformation for the twinning-buckling of tatb under dynamic loading[END_REF][START_REF] Lafourcade | Elastic instability in graphite single crystal under dynamic triaxial compression: Effect of strain-rate on the resulting microstructure[END_REF][START_REF] Leyssale | Mechanisms of elastic softening in highly anisotropic carbons under in-plane compression/indentation[END_REF]. Shear strains are applied in the [-1.0, +1.0] % range with the same increment. After the strain is applied, the simulation cells are equilibrated in the NVT ensemble for 1.0 ns, leading to the stress relaxation of the the system through the relaxation of atoms positions. Stress tensors are then obtained by averaging stress components over the last 100 ps of the trajectories and isothermal elastic constants C T ij are determined via linear fits of the corresponding stress vs. strain curves.

Adiabatic elastic constants

Adiabatic elastic constants are obtained in a somehow opposite way. First, the system is equilibrated for 1 ns in the NVT ensemble using the equilibrium lattice parameters obtained from the NPT trajectory at 1 atm. Then, 20 independent configurations are selected from the last 100 ps of the NVT simulation, on which finite adiabatic strains, spanning the same values as those used in the isothermal case, are applied via pure homothethic transformations to the simulation cells. The resulting 20 uncorrelated stress-strain curves, for each deformation type, are then gathered to build the average stress-strain curve, from which the corresponding adiabatic, or isentropic elastic components C S ij are obtained, where the superscript S relates the elastic constant to entropy. This way, one deformation type gives directly access to 6 elastic constants.

Results

Room temperature lattice parameters

The lattice parameters of graphite single crystal at ambient temperature, i.e. P = 1 atm and T = 300 K, predicted using the different potentials are given in

Room temperature elastic properties

Isothermal elastic constants computed at P = 1 atm and T = 300 K using the different potentials are compared in Table 2 to experimental values from two reference papers. It has been previously presented that some potentials may slightly deviate from the expected transverse isotropy symmetry of graphite when regarding elasticity [START_REF] Farbos | Nanoscale elasticity of highly anisotropic pyrocarbons[END_REF][START_REF] Gamboa | On the prediction of graphene's elastic properties with reactive empirical bond order potentials[END_REF]. Hence, as discussed before, we adopt in the present work the general case of orthotropic symmetry of the second-order elastic tensor in order to compare the computed elastic constants of graphite single crystal.

Table 2: Isothermal elastic constants (in GPa) of graphite single crystal at P = 1 atm and T = 300 K computed using the different interatomic potentials under the assumption of an orthotropic elastic symmetry. The last column corresponds to the metric d in ti that is meant to measure the deviation from the transverse isotropy, computed from the in-plane components of the second-order elastic tensor. However, more recent theoretical [START_REF] Cousins | Elasticity of carbon allotropes. iii. hexagonal graphite: Review of data, previous calculations, and a fit to a modified anharmonic keating model[END_REF] and experimental investigations [START_REF] Grimsditch | Shear elastic modulus of graphite[END_REF] suggest that the recommended by Cousins and Heggie [START_REF] Cousins | Elasticity of carbon allotropes. iii. hexagonal graphite: Review of data, previous calculations, and a fit to a modified anharmonic keating model[END_REF].

Two metrics are defined to identify the deviation from transverse isotropy symmetry, regarding in-plane and our-of-plane symmetry, respectively:

d in ti = 100 ×          1 - 1 2 ×          C 2 11 + C 2 22 2C 11 C 22 + 2C 66 2 C 2 11 +C 2 22 C 11 +C 22 -C 12                   (3)
d out ti = 100 × 1 - 1 2 × C 2 13 + C 2 23 2C 13 C 23 + C 2 44 + C 2 55 2C 44 C 55 (4)
where each ratio should be equal to 1 in the presence of an elastic tensor with transverse isotropy symmetry.

The values of d in ti are reported in the last column of Table 2, noticing that d in ti = 0 is expected experimentally. d out ti is not reported for the isothermal constants due to the very low values of C 44 , and the necessarily large associated relative uncertainty. Amongst the different potentials LCBOPII is the one that produces the second-order elastic tensor with the closest symmetry to in-plane transverse isotropy, followed by CEDIP, the three other potentials showing much larger deviations.

In Table 3 are listed the adiabatic elastic constants computed with the different potentials. Since no experimental data correspond to adiabatic elastic constants measurements for graphite single crystal, we only discuss them in terms of their differences with isothermal elastic constants as well as their deviation from the transverse isotropy symmetry. As a fectly identical to the the experimental data of Blakslee et al. [START_REF] Blakslee | Elastic constants of compressionannealed pyrolytic graphite[END_REF], which indicates that these potentials have been calibrated to reproduced the elastic constants of graphite using homogeneous deformation without relaxation of atomic positions, i.e. under adiabatic conditions. In addition, all the potentials lead to a value of the adiabatic constant C 66 that is higher than the isothermal constant, whereas the adiabatic constant C 12 is lower than the isothermal value. Under the assumption that graphite single crystal is transversely isotropic, the relationship C 66 = (C 11 -C 12 )/2 holds and a decrease in C 12 systematically involves an increase in C 66 , which could explain the observed trend.

Temperature evolution of lattice parameters

Figure 1 compares the evolution with temperature of the lattice parameters computed using the different potentials to available experimental data [START_REF] Nelson | The thermal expansion of graphite from 15 c. to 800 c.: part i. experimental[END_REF][START_REF] Kellett | A study of the amplitude of vibration of carbon atoms in the graphite structure[END_REF][START_REF] Walker | X-ray diffraction studies of a graphitized carbon -changes in interlayer spacing and binding energy with temperature[END_REF][START_REF] Kellett | The c-axis thermal expansion of carbons and graphites[END_REF]. Common to all the potentials is the observation of a monotonous increase in the in-plane lattice parameter (a) with temperature (see Figure 1a), which conflicts with experimental results

showing that a passes through a minimum at about 700 K. The slight contraction observed experimentally at low temperatures can indeed not be captured by classical MD as pointed out in recent investigations comparing classical and quantum, path integral, MD in the case of monolayer graphene [START_REF] Brito | Quantum effects in a free-standing graphene lattice: Path-integral against classical monte carlo simulations[END_REF][START_REF] Herrero | Quantum effects in graphene monolayers: Path-integral simulations[END_REF]. Besides this, all potentials seems to significantly overestimate both in-plane and out-of-plane thermal expansion, even at large temperatures, CEDIP being the closest to experimental data. Quantitatively, predicted a values using REBO, CEDIP and LCBOPII are very close to the commonly accepted experimental data of Kellet and Richards [START_REF] Kellett | The thermal expansion of graphite within the layer planes[END_REF] at low temperatures and differ by less than ∼ 1 % at T = 3000 K. Conversely, the values obtained for AIREBO and AIREBO-M are significantly lower, as already discussed above for the room temperature data. Due to the non-linear evolution of graphite single crystal elasticity (both isothermal and adiabatic) with temperature, we propose an analytical expression for this dependency based on polynomial fits. For the sake of simplicity, only the LCBOPII potential is considered in the remaining of the manuscript. Indeed, this potential was shown to produce reasonable lattice parameters and elastic constants, preserving well the trnsverse isotropy of graphite.

Furthermore, it allows to capture out-of-plane elastic constants, which CEDIP and REBO cannot. We take the condition T r ef = 300 K as the reference configuration, leading to the following formulation of the isothermal/adiabatic spatial elasticity tensor along the ambient pressure isobaric pathway:

c S/T (T ) = c S/T r ef (T r ef ) I + N k=1 A k (T -T r ef ) k (5)
where c S/T r ef is the stiffness tensor at the reference state, with superscripts S and T corresponding to isothermal and adiabatic forms, respectively. A k is a tensor containing the polynomial fit coefficients and N is taken equal to 3. It is to be noted that the product between c S/T r ef and A k is done through the Hadamard product, i.e. component-wise.

The temperature-dependent isothermal and adiabatic elastic coefficients and their corresponding polynomial fits are represented in Figure 4. In addition, the polynomial coefficients for both isothermal and adiabatic elastic tensors are provided in the supplementary material.

Temperature dependence of elastic anisotropy

When dealing with highly anisotropic materials like graphite, different metrics computed from the elastic tensor can help to discriminate whether the proposed analytical formulation is reliable or not. One of these metrics is the universal anisotropy index [START_REF] Ranganathan | Universal elastic anisotropy index[END_REF] which computes the level of anisotropy of any material from its elastic constants. It is based on the bulk and shear moduli computed using the Voigt and Reuss approximation. The Voigt bound relates to the behavior of a material under an isotropic deformation-imposed process while the Reuss bound is related to the behavior of a material under an imposed hydrostatic pressure. It is obvious that any isotropic material will behave in the same way under both conditions. Yet, this is not the case in the presence of anisotropy. The universal anisotropy index reads:

A U = c : s = 5 G V G R + K V K R -6, (6)
which takes a value of 0 for an isotropic material, since Voigt and Reuss estimates for bulk (K) and shear (G) modulus are identical. Figure 5 shows the evolution with temperature of A univ along the 1 atm isobaric pathway computed from MD data and the isothermal and adiabatic analytical expressions.

First, we see that analytical formulations reproduce well both isothermal and adiabatic MD data, thus validating the polynomial fits. More importantly, Figure 5 shows that the elastic anisotropy is severely temperature dependent. At ambient temperature, anisotropy is very high compared to a large panel of materials [START_REF] Ranganathan | Universal elastic anisotropy index[END_REF][START_REF] Kube | Elastic anisotropy of crystals[END_REF] and takes values of approximately 93 and 2169 in adiabatic and isothermal conditions, respectively. Concerning the evolution with temperature, the isothermal anisotropy decreases non-linearly with temperature and is divided by almost a factor of 5.40 at 4000 K with respect to the 300 K case. The adiabatic anisotropy follows a very similar evolution and the ratio between the 300 K and the 4000 K values is around 4.65.

3.6. Adiabatic vs. isothermal elastic constants: relationship and implications We show in Figure 6 the ratios of the isothermal over adiabatic elastic constants as a function of temperature. Interestingly, these ratios remain almost constant for C 11 , C 66 and C 44 , with values around 0.8 for C 11 and C 66 , and about 0.02 for C 44 , although a slight linear decrease with increasing T is observed, especially for C 11 . Here it is important to note that any value different from unity of the

C T ij C S ij
ratio indicates that stress relaxation will take place during, or after, deformation, depending on the deformation rate. For instance, the value of ∼ 0.8 obtained for C 11 indicates that when single crystalline graphite is subjected to a fast uniaxial in-plane deformation, the increase in stress is systematically followed by a relaxation process. In other words, this means that the measured elastic tensor will depend on strain rate, with the adiabatic and isothermal constants constituting the instantaneous and quasistatic limits, respectively. Similarly, the extremely low value of the isothermal over adiabatic ratio for C 44 implies a systematic accommodation of interlayer shear deformation at large times.

Conversely to the formerly discussed constants, the isothermal over adiabatic ratio significantly decreases with T, with a close to linear dependency, for the three other constants, indicating that stress relaxation is strongly temperature-dependent for these constants. Ratios for C 12 and C 13 both decrease from about 1.5 at 300 K to about 0.8 at 4000 K, while for C 33 a decrease from about 0.85 at 300 K to about 0.15 at 4000 K is observed. In order to understand the transition from adiabatic to isothermal elastic response, we used MD simulations with the LCBOPII potential to study the mechanical response of a graphite single crystal at various temperatures and strain amplitudes. A graphite single crystal of approximate dimensions 10 × 10 × 10 nm 3 containing 107520 atoms was equilibrated in the NVT ensemble for 50 ps for temperatures between 300 K and 4000 K. After the material is fully equilibrated, it is instantaneously dilated along the [100] direction for different strains ε between 0.1 % and 1 % in order to remain within the elastic domain, similar to the one used to compute the elastic constants. This deformation process is performed without any thermostat coupling in order to let the material relax over time without impacting the relaxation kinetics, which might be related to the thermostat damping constant. However, the deformation induced relaxation systematically involves heating of the system due to the total energy conservation. Yet, this process allows capturing the dynamics of the elastic response. Indeed, when the deformation is applied, the stress level is directly related to the adiabatic elastic constant C S 11 at the temperature right before the deformation, while after the relaxation process, the stress level corresponds to the isothermal elastic constant C T 11 at the temperature at the end of the relaxation process. What happens in between is then directly related to the transition between adiabatic and isothermal elastic response. Performing this procedure at different initial temperatures and different strains allows to identify how these conditions affect the transition, which can be related to the relaxation process behind the front of a shock-wave for example. In Figure 7, we display the time evolution of the stress component σ 11 , for different longitudinal strains ε 11 up to 1 %. The different blue and cyan dashed lines correspond to the analytical limits σ ∞ 11 and σ 0 11 computed from isothermal and adiabatic elastic constants using linear elasticity:

σ 0 11 = C S 11 (T 0 ) • ε 11 , (7)
σ ∞ 11 = C T 11 (T ∞ ) • ε 11 . (8)
where T 0 and T ∞ correspond to the temperatures right before the deformation is applied and at the end of the relaxation process, respectively. A focus on the case ε 11 = 0.01 is displayed in Figure 7, where the good agreement between the stress levels and the analytical values can be seen. However, one can see that the stress is oscillating around σ ∞ 11 while decreasing in amplitude over time. This oscillation of the stress is due to the elastic energy stored in the material at the instantaneous deformation that is relaxing through elastic waves emissions.

In order to extract a characteristic time from this temporal evolution, the first step consists in extracting the envelope of the stress signal, that is directly related to the decrease of the stress over time. We then consider the following exponential decay function that relates the stress versus time to the isothermal and adiabatic elastic constants using the previously defined analytical limits:

σ 11 (t) = σ ∞ 11 + (σ 0 11 -σ ∞ 11 ) exp(-t/τ ) (9)
with τ the characteristic time associated to the exponential decay. We then define the characteristic relaxation time t r as the time at which the decrease is completed at 95 %.

The relaxation time for an initial temperature of 300 K and ε = 0.01 is displayed as a vertical gray line on Figure 7b. This procedure has been applied to stress signals for different strains ε 11 ∈ [0.001, 0.1] and for temperatures T ∈ [300 K, 4000 K]. Figure 8a shows that the relaxation time t r for instantaneous in-plane tensile loads is in the ps range for the considered strain and temperature range and that it increases with applied strain while it decreases with temperature. In addition, it is also interesting to discuss elastic relaxation in terms of its characteristic length:

L r (T, ε) = t r (T, ε) C T 11 (T )/ρ(T ) (10)
where the square root term corresponds to the longitudinal elastic wave velocity of graphite, computed from the ratio between the isothermal elastic constant C T 11 (T ) and the density ρ(T ) at finite temperature. As shown in Figure 8b, similarly to t r , L r evolves slightly with temperature, reaching sort of a plateau at large strains. Indeed, the lower bound of this characteristic length remains below 10 nm at ε 11 = 0.01 for temperatures above 1000 K while it is always below 25 nm for the entire set of simulations. This implies that, when a tensile (for example after a reflection) shock-wave propagates along [100] direction in a graphite single crystal, a slab of width < L r = 15 nm behind the front is required to accommodate the local high strain-rate and to transit from adiabatic to isothermal conditions. When considering high strain-rate deformation, one can talk in terms of relaxation time instead of characteristic length, since no discontinuities of particle velocities exist in the material. This would mean that if the dynamic deformation process involves local deformations above 1 % within a time window lower than a few ps, the process could be considered as adiabatic. On the other hand, a transition from adiabatic to isothermal would be involved.

Discussion and conclusion

Equilibrium MD simulations of graphite single crystal have been performed on an extended thermodynamic domain (up to 4000 K) at ambient pressure using various potentials of the literature, namely REBO, AIREBO, AIREBO-M, CEDIP and LCBOPII. The lattice parameters evolution with temperature has been computed with the different models and the LCBOPII model has been identified as a good candidate to represent the thermal behavior of graphite single crystal. Additionally, at each temperature, the spatial isothermal and adiabatic elasticity tensors have been computed through the derivation of the Cauchy stress with respect to strain, using the different models considered in this work. Adiabatic elastic constants correspond to the instantaneous response of a material to a mechanical deformation of very short wavelength, i.e. as found in the vicinity of the shock-wave front or during very high strain-rate deformation. On the opposite, isothermal elastic constants correspond to a fully equilibrated mechanical response at longer time scales, i.e. after a shock-wave has passed or under quasi-static or low strain-rate deformations.

Through the entire temperature range, both isothermal and adiabatic elastic coefficients are found to evolve non-linearly with temperature, independently from the potential used.

Based on the good agreement with the literature for both lattice parameters evolution with temperature and for ambient temperature isothermal elastic constants, the LCBOPII model seems appropriate to study the mechanical behavior of graphite at high temperature. An analytical formulation for both isothermal and adiabatic spatial elasticity tensors as a function of temperature has been introduced as a preliminary element in order to build a mesoscopic model for the elastic behavior of graphite single crystal. The dependence on temperature is reproduced through a component-wise third-order polynomial fit using the Hadamard product, allowing to interpolate the elastic coefficients calculated using MD simulations.

Interestingly, it has been observed that isothermal and adiabatic elastic constants significantly differ from each other. A clear consequence of that is the difference in the elastic anisotropy of graphite, the universal anisotropy index being larger by a factor of ∼ 20 in the isothermal case than in the adiabatic case, whatever the temperature. Also, temperature has a tremendous influence on elastic anisotropy with a decrease by a factor of five of the universal anisotropy index observed between 300 K and 4000 K. This is of high importance since it has been shown that elastic softening deformation mechanisms such as elastic instabilities and in particular buckling are driven by the material anisotropy [START_REF] Lafourcade | Irreversible deformation mechanisms for 1,3,5-triamino-2,4,6trinitrobenzene single crystal through molecular dynamics simulations[END_REF][START_REF] Lafourcade | Mesoscopic constitutive law with nonlinear elasticity and phase transformation for the twinning-buckling of tatb under dynamic loading[END_REF][START_REF] Lafourcade | Elastic instability in graphite single crystal under dynamic triaxial compression: Effect of strain-rate on the resulting microstructure[END_REF][START_REF] Leyssale | Mechanisms of elastic softening in highly anisotropic carbons under in-plane compression/indentation[END_REF]. The presented results actually suggest that such buckling transition may be suppressed at high temperatures and/or under adiabatic conditions.

Finally, MD simulations were used to investigate the details of the adiabatic to isothermal elastic relaxation by imposing an instantaneous deformation to a graphite single crystal at different temperatures and strain levels. The simulations suggest that in the presence of a tensile shock-wave going through a graphite crystal, a non negligible width of ∼ 20 nm is required to accommodate a 1 % instantaneous in-plane tension applied to graphite at room temperature and to allow the local material state to transit from adiabatic to isothermal stress response. Overall, this study provided a detailed, temperature dependent, analytic model of the adiabatic and isothermal elasticity of single crystalline graphite, that can be integrated into mesoscale models of graphite-based materials elasticity.

value of 5

 5 GPa is correct for C 44 , the lower values obtained on compressed pyrolytic graphite by Blakslee et al. being attributed to glissile basal dislocations. Therefore, it seems that all the considered potentials underestimate C 44 . As for C 13 , the three potentials predict values in the 6-8 GPa range, that are intermediates to the values of Blakslee et al.[START_REF] Blakslee | Elastic constants of compressionannealed pyrolytic graphite[END_REF] and Bosak et al.[START_REF] Bosak | Elasticity of single-crystalline graphite: Inelastic x-ray scattering study[END_REF], 15 and 0 GPa, respectively, and almost perfectly equal to the value of 7.9±3.5

Figure 1 : 3 . 4 .Figure 2 2 , C 13 +C 23 2 and C 44 +C 55 2 ,

 13422132442 Figure 1: Evolution of (a) in-plane and (b) out-of plane lattice parameters with temperature obtained computed with the different potentials. Experimental data from Nelson and Riley [56], Kellet and Richards [62, 59], Walker, McKinstry and Wright [58] and Kellett, Jackets and Richards [57] are given for comparison.

Figure 2 :

 2 Figure 2: Evolution with temperature of the isothermal elastic constants of graphite single crystal at P = 1 atm computed using the different potentials. For clarity, error bars are only given for data obtained with the LCBOPII potential. Error bars for the other data sets are expected to be of the same order. Dashed lines for C 66 data correspond to predictions obtained from C 11 and C 12 assuming transverse isotropy symmetry.

Figure 3 :

 3 Figure 3: Same as Figure 2 for the adiabatic elastic constants.

Figure 4 :

 4 Figure 4: Comparison of MD data and polynomial fits for the temperature dependence of the isothermal and adiabatic elastic constants. symbols: MD data; lines: polynomial fits.

Figure 5 :

 5 Figure 5: Universal Anisotropy Index as a function of temperature computed from (a) isothermal and (b) adiabatic elastic tensors, comparison between analytical formulation and MD data. All data were obtained with the LCBOPII potential.

Figure 6 :

 6 Figure 6: Evolution of the ratio between isothermal (C ij T ) and adiabatic (C S ij ) elastic constants of graphite single crystal using the LCBOPII potential as a function of temperature. All data were obtained with the LCBOPII potential.

Figure 7 :

 7 Figure 7: Evolution of (a) the stress component σ 11 vs time after an instantaneous tensile strain ε 11 is applied to the sample, for different values of ε 11 . The blue and cyan dashed lines represent the limit values, corresponding to the analytic stress level using the isothermal and adiabatic C 11 coefficients, respectively. (b) Focus on the case ε 11 = 0.01. The dark red line is the computed limit envelope of the stress signal over time while the pink dashed line is the result of the fitted exponential decay function, with characteristic time τ . The vertical gray line corresponds to the relaxation time t r computed as described in the text.

Figure 8 :

 8 Figure 8: Evolution with strain and temperature of the characteristic (a) time and (b) length of the adiabatic to isothermal elastic stress relaxation after uniaxial in-plane compression. All data were obtained using the LCBOPII potential.

Table 1

 1

	. Experimental

AIREBO and AIREBO-M give a better prediction of the inter-layer spacing (d 002 = 0.5c) while LCBOPII overestimates it by 1.6 %. As discussed above, no c values could be predicted using REBO and CEDIP.

Table 1 :

 1 Lattice parameters of graphite at P = 1 atm and T = 3000 K computed using the different potentials and compared to the experimental reference[START_REF] Baskin | Lattice constants of graphite at low temperatures[END_REF]. Relative errors in % with respect to the reference are reported in parentheses.

		a (Å )	c (Å )
	Exp.	2.4590	6.7076
	LCBOPII	2.4590 (0.000) 6.8180 (1.646)
	AIREBO	2.4193 (-1.615) 6.7078 (0.004)
	AIREBO-M 2.4193 (-1.616) 6.7078 (0.003)
	REBO	2.4600 (0.042) N/A
	CEDIP	2.4605 (0.059) N/A

  In a general way, all the force fields respect relatively well the hierarchy of the different elastic constants with C 11 > C 66 > C 12 > C 33 > C 13 > C 44 . The potentials also perform reasonably well from the quantitative perspective. Albeit being slightly lower than the experimental references, the error on computed C 11 ranges from about 10 % with AIREBO to about 27 % with REBO. As for C 11 , C 66 is underestimated by 8 to 15 % with AIREBO and AIREBO-M, the largest error, 24-30 %, being obtained for REBO and CEDIP. Regarding C 12 the value predicted with LCBOPII, 159.9 GPa, is almost in the middle of the two experimental reference values, the other potentials predicting slightly larger values, the largest, 253 GPa, obtained with CEDIP being 40-80 % too large. Predicted C 33 values are either about 7-8 GPa larger (AIREBO) or lower (LCBOPII, AIREBO-M) than the experimental data.Concerning the shear elastic constants C 44 and C 55 , there is a large experimental spreading of the data since Bosak et al.[START_REF] Bosak | Elasticity of single-crystalline graphite: Inelastic x-ray scattering study[END_REF] and Blakslee et al.[START_REF] Blakslee | Elastic constants of compressionannealed pyrolytic graphite[END_REF] predicted values of 5 GPa and 0.18 GPa, respectively. LCBOPII, AIREBO and AIREBO-M potentials lead to consistent values in the 0-0.2 GPa range, i.e. close to the experimental data of Blakslee et al.[START_REF] Blakslee | Elastic constants of compressionannealed pyrolytic graphite[END_REF].

		C 11	C 22	C 33	C 12	C 13	C 23	C 44	C 55	C 66	d in ti
	Exp. [21]	1060.0	36.5 180.0	15.0	0.18	440.0	0.0
	Exp. [22]	1109.0	38.7 139.0	0.0		5.0		485.0	0.0
	LCBOPII	941.0	945.0 28.5 159.9 6.5	6.1 0.23 0.17 389.2 0.05
	AIREBO	1000.9 1000.4 45.6 209.9 7.6	6.5 0.15 0.05 407.8 -1.54
	AIREBO-M	936.5	913.1 29.6 190.4 7.6	7.9 0.04 0.01 406.2 -3.60
	REBO	808.7	827.7 N/A 183.2 N/A N/A N/A N/A 340.0 -5.19
	CEDIP	920.6 921.40 N/A 252.9 N/A N/A N/A N/A 334.6 -0.15

Table 3 :

 3 Adiabatic elastic constants (in GPa) of graphite single crystal at P = 1 atm and T = 300 K computed using the different interatomic potentials under the assumption of an orthotropic elastic symmetry. The last two columns correspond the in-plane d in ti and out-of-plane d out and for all the potentials, the computed adiabatic elastic constants are much closer to the transverse isotropy symmetry than their isothermal counterparts, with one noticeable exception being the large C 11 /C 22 ratio, 1.146, obtained with REBO. One notice that the adiabatic values of C 11 obtained with LCBOPII and REBO are almost per-

	ti	metrics measuring the deviation from transverse

with 5 independent elastic component. In the following, we consider the general case with elastic orthotropic symmetry, in order to compare the elastic constants computed with different potentials, and determine whether they verify the transverse isotropic symmetry. By construction, it is necessary to apply 6 deformation types to get the full second-order elastic tensor from MD simulations. As for the lattice constant calculations, elastic constants involving the out-of-plane direction are not investigated with REBO and CEDIP. Mechanical properties are computed by applying an homothetic strain to the entire system with atomic positions being remapped into the box, followed or not by an equilibration in the NVT ensemble whether we aim at calculating isothermal or adiabatic elastic constants. The differences between the two methods are shortly explained in the following sections.

Acknowledgments

The authors would like to thank Dr. Christophe Denoual and Dr. Nicolas Pineau for their help in both setting up the MD simulations using the LCBOPII potential and discussing the results. This research is financially supported by the French Ministry of Defence -Defence Innovation Agency. Atomistic simulations were performed using the computational resources of the TGCC-CEA.


Tingting Zhai

Matteo Mazzola

Demetrio Macias

Rafael Salas-Montiel
email: rafael.salas@utt.fr

Efficient spin-direction coupling between circular polarized electric and magnetic dipoles and photonic modes

Integrated single photon sources are essential elements in quantum technology applications such as computing, simulation, communication, metrology and in quantum photonic neural networks. The directional radiation and scattering of integrated single-photon sources play a crucial role in light manipulation, which relies on the electric and magnetic dipole moments of the emitter sources. Although clever physical insights and designer intuition strategies have been successfully applied in the development of integrated sources, inverse design strategies could enhance its performance. Inverse strategies such as topology optimization produce optimized geometries for nanophotonic devices that maximize its performance. Recently, topology-optimized couplers for on-chip single-photon sources were designed to efficiently couple light between a photonic mode in a waveguide and a single-photon source, assumed to be an electric dipole oriented along one axis. However, the superposition of orthogonal electric and magnetic dipoles can also be harnessed due to the additional degrees of freedom via their interference. In this work, we extended the topology optimization to couplers that enhance spin-direction coupling of circularly polarized, Huygens, and Janus dipoles. We demonstrated that topology optimization not only increases the light coupling to a desire mode but also enhances the electric and magnetic local density of states by the modification of the electromagnetic environment while maintaining the amplitude and phase relation between the orthogonal dipoles for unidirectional coupling. Currently, coupling efficiency and enhanced LDOS of up to 88%, 94%, and 93 % and up to 9.4, 18.6, and 14.9 are obtained for the circular polarized, Huygens, and Janus dipoles, respectively. The adjoint topology optimization improves the performances of 3D integrated photonic devices.

Introduction

Integrated single photon sources are essential integrated photonic components within a full integrated quantum photonic circuit. They provide the necessary generation of individual photons for further single and multiple states processing and detection. Such components will form the next generation of photonic circuits that will open the way to quantum photonic technologies. Practical quantum photonic technologies such as computing, simulation, communication, and metrology as well as recent advances in quantum optical neural networks are envisioned with the use of integrated photonic circuits [START_REF] Steinbrecher | Quantum optical neural networks[END_REF]. Among the single photon sources available, quantum emitters are promising sources for quantum networks [START_REF] Lu | Quantum-dot single-photon sources for the quantum internet[END_REF] as they can reach near-optimal emission production of on-demand and single indistinguishable photons [START_REF] Somaschi | Near-optimal single-photon sources in the solid state[END_REF]. In addition to quantum dot emitters, quantum emitters embedded in 2D transition metal dichalcogenides have shown to present unique properties [START_REF] Kianinia | Quantum emitters in 2D materials: Emitter engineering, photophysics, and integration in photonic nanostructures[END_REF] and have been integrated in photonic structure [START_REF] Peyskens | Integration of single photon emitters in 2D layered materials with a silicon nitride photonic chip[END_REF]. It is known that the rate of emission of photons by quantum emitters is not an intrinsic property and can therefore be enhanced via the modification of their local environment or equivalently, via the modification of the local density of states (LDOS) [START_REF] Purcell | Spontaneous emission probabilities at radio frequencies[END_REF]. The coupling of quantum emitters to photonic microcavities with optimized spatial and spectral overlap has been proven to increase the rate of emission [START_REF] Fröch | Coupling hexagonal boron nitride quantum emitters to photonic crystal cavities[END_REF]. In addition to the coupling of quantum emitters to photonic structures that enhance the Purcell effect, the emission radiation is directed into a particular direction for further processing [START_REF] Arcari | Near-unity coupling efficiency of a quantum emitter to a photonic crystal waveguide[END_REF][START_REF] Olthaus | Optimal photonic crystal cavities for coupling nanoemitters to photonic integrated circuits[END_REF]. In this regard, integrated photonics is an excellent platform for both increasing the LDOS and collection efficiency. Quantum dots emitters coupled to planar nanobeam waveguide [START_REF] Kirsanske | Indistinguishable and efficient single photons from a quantum dot in a planar nanobeam waveguide[END_REF] and electrically contacted planar photonic crystal waveguide membrane has demonstrated to produce indistinguishable and efficient single photons [START_REF] Uppu | Scalable integrated single-photon source[END_REF]. Hybrid structure consisting of a high-index layer placed on top of an ion-exchanged glass waveguide has also proved the enhancement of emission of quantum emitters [START_REF] Madrigal | Hybrid integrated optical waveguides in glass for enhanced visible photoluminescence of nanoemitters[END_REF]. Moreover, plasmonic waveguides have also been employed for the enhancement of the emission of quantum emitters in quantum technologies [START_REF] Kumar | Single photon emitters coupled to plasmonic waveguides: A review[END_REF]. However, all these integrated photonic structures were developed following physical insights and designer's intuition, dealing with predefined configurations. As the quantum integrated circuits scale up to large bandwidth, multi-frequency applications, nonlinear phenomena, and dense integration, the prototypical approach poses a challenge of increasing complexity [START_REF] Molesky | Inverse design in nanophotonics[END_REF]. Therefore, controlling the propagation of guided light via fast and robust near field interference between polarization components of a single photon sources is required [START_REF] Picardi | Experimental demonstration of linear and spinning Janus dipoles for polarisation-and wavelength-selective near-field coupling[END_REF]. As an alternative candidate, inverse design strategies such as topology optimization could reshape the landscape of available nanophotonic structures that have more forgiving design rules with the goal of maximizing the performance [START_REF] Jensen | Topology optimization for nano-photonics[END_REF]. Recently, a nanoscale light extractor was designed using adjoint-optimization methods to enhance by 35 times the broadband emission of nitrogen vacancy (NV) centers in diamond. The extractor collects the emission into radiative modes (i.e., to free space), directing the light into a narrow cone that can be collected with low numerical aperture optical systems with 40% collection efficiency [START_REF] Chakravarthi | Inverse-designed photon extractors for optically addressable defect qubits[END_REF]. Another extractor, designed using inverse design, enhances by 14 times the emission of NV centers in diamond [START_REF] Wambold | Adjoint-optimized nanoscale light extractor for nitrogen-vacancy centers in diamond[END_REF]. Topology-optimized couplers for on-chip single-photon sources were designed to efficiently couple light emission of a single-photon emitter to a photonic mode in a waveguide. The emitter is however, assumed to be an electric dipole oriented along one axis [START_REF] Yesilyurt | Efficient topologyoptimized couplers for on-chip single-photon sources[END_REF].

Circular polarized dipoles, together with Huygens and Janus dipoles, form the complete set of directional dipolar sources in the far-and near-field [START_REF] Picardi | Janus and Huygens dipoles: Near-field directionality beyond spin-momentum locking[END_REF]. These dipoles are known to have a near-field directionality and can thus excite a guided mode in one direction only if placed in the vicinity of a waveguide. The near field directionality for the circularly polarized dipole is explained based on its spin angular momentum. This dipole is formed by the superposition of two orthogonal electric or magnetic dipole, 90° out of phase. The transverse spin is at the origin of the directional excitation of guided modes via spin-orbit coupling. This coupling is done by placing the dipole in regions where the modes show local circular polarization [START_REF] Espinosa-Soria | Transverse Spin and Spin-Orbit Coupling in Silicon Waveguides[END_REF]. In the Huygens dipole, an electric and a magnetic dipoles are orthogonal to each other and in phase while in the Janus dipole the orthogonal electric and magnetic dipoles have a 90° phase shift. These combinations provide directionality in both the far-and near-fields. The Janus dipole is said to have two "faces": one face couples into the mode, while the other is noncoupling [START_REF] Michela | Integrated Janus dipole source for selective coupling to silicon waveguide networks[END_REF]. Although these dipoles are highly directional and have unique properties for side coupling to guided modes, their coupling efficiency is low, and their local densities of states (LDOS) are weakly enhanced due to the side coupling.

In this work, we thus extended the topology optimization to couplers with circularly polarized dipoles as well as Huygens and Janus dipoles. We demonstrated that topology optimization not only increases the light coupling to a desire mode of an output waveguide but also it enhances the LDOS. Specifically, we calculate the coupling efficiency and LDOS enhancement of a topology optimized structure compared to a non-optimized structure.

Methods and Materials

We used the adjoint method implemented by Meep [START_REF] Oskooi | Meep: A flexible free-software package for electromagnetic simulations by the FDTD method[END_REF] to systematically design a compact silicon photonic structure that maximize a fraction of the radiation emission of a single-photon emitter into the fundamental mode of the output silicon waveguide (Fig. 1). The single-photon emitter is modeled as a dipole as described in Table I and thus, we designed the structure to maximize the coupling efficiency into the fundamental mode of the output waveguide.

Table I. Used sources and their respective electric and magnetic dipoles.

Circular polarized source

Huygens' source

Janus source noncoupling coupling Electric dipole p = (1, 0, 1i) p = (0, 0, 1)

p = (-1, 0, 0) p = (+1, 0, 0) Magnetic dipole m = (0, 0, 0) m = (0, 1c, 0) m = (0, 1ic, 0) m = (0, 1ic , 0)
In the simulations, the computational domain is a rectangle of size 13 µm by 7 µm, including 1-µm-thick perfectly matched layers (PML) around the four boundaries. We use a grid mesh of 40 (pixels.μm -1 ), that is 1/40=0.025 µm. The optimization domain is 3 µm x 3 µm and has thus 120 by 120 design parameters. The right-hand side of the design region is connected to a photonic waveguide of width 0.5 µm and length 4 µm and refractive index n2 =3.4. The outside region is silicon dioxide with index of refraction n1 = 1.44 [START_REF] Kischkat | Mid-infrared optical properties of thin films of aluminum oxide, titanium dioxide, silicon dioxide, aluminum nitride, and silicon nitride[END_REF] (Fig. 1b). We specify a characteristic length scale of the structure in the domain a = 1 µm. Then, a wavelength of λ = 1.55 μm corresponds to a frequency 0.64516 (i.e., a/λ = 1μm/1.55μm = 0.64516). Initially a dipole is placed at the center of the design region.

The objective function is defined by:

𝜂 𝑖 (𝐫 𝟎 , 𝜔) = |𝑎 𝑖 (𝐫 𝟎 ,𝜔)| 2
𝑃 𝑙,𝑖𝑛 (𝐫 𝟎 ,𝜔)

…. [START_REF] Steinbrecher | Quantum optical neural networks[END_REF] ,

where ai(r0, ω) is the expansion coefficient of the i-th mode at the output waveguide. The expansion coefficient is given by an overlap integral over a cross-section A, ai(r0, ω) = < Ωi | Ω >, where Ωi is the eigenmode profile corresponding to the i-th eigenmode at the output waveguide and Ω is the total fields. Pl,in is the power emitted by the dipole placed at r0 and oriented along the direction 𝑙 ∈ {x, y, z}. This input power was calculated by calculating the Poynting flux power transmitted out of a box surrounding the dipole.

The optimization is done within a spectrum centered at frequency 0.64516 with a frequency width of 0.129, from 0.58 to 0.71 in units of (a/λ). This means at the central wavelength λ = 1.55 μm, from 1.408 μm to 1.724 μm. In other words, we seek to maximize the fraction of the power emitted by the source into the fundamental mode of the output waveguide in a spectrum centered at 0.64516 with a frequency width of 0.129 in units of (a/λ).

In addition, we calculated the modification of the LDOS of the source placed in the topology-optimized structure. The total LDOS(r0, ω) is given by the sum of the electric and magnetic contribution [START_REF] Joulain | Definition and measurement of the local density of electromagnetic states close to an interface[END_REF]. In FDTD, the electric LDOS is given by [START_REF] Oskooi | Advances in FDTD computational electrodynamics: photonics and nanotechnology[END_REF]:

LDOS l (𝐫 𝟎 , 𝜔) = - 2 𝜋 𝜖(𝐫 𝟎 )𝑅𝑒 [ 𝐸 𝑙 (𝐫 𝟎 ,𝜔)𝐩(𝜔) * |𝐩(𝜔)| 2 ], … (2)
where the subscript "𝑙" indicates the oriented of the dipole, 𝑙 ∈ {x, y, z} and |p(ω)| 2 is the spectrum of an impulsive point-dipole source 𝐉(𝐫 𝟎 , 𝑡) = 𝐞 𝐥 𝛿(𝐫 -𝐫 𝟎 )𝑝(𝑡). However, this value in the denominator is cancelled as we are not interested in the absolute value of the LDOS, but rather in the relative enhancement of the LDOS relative to another structure. We thus obtained a normalized LDOS by dividing the LDOS in the optimized region by the LDOS of the same dipole placed in a non-optimized region (i.e., homogeneous design region with index n2 = 3.4). This ratio is equivalent to the Purcell factor.

The size of the design region of 3 µm x 3 µm is chosen so that the dipole has multiple modes to interact with. The dipole interacts with the non-zero field components ez, hy, and ex of the available modes to perform directionality and optimization. The optimization algorithm forces the amplitude and phase relationship of the dipoles to maximize the coupling efficiency into the fundamental mode. In our case, the fundamental mode has the main component of the magnetic field out-of-plane (i.e., component Hy).

The dispersion relation of the output waveguide is plotted in Fig. 1d. At the central frequency 0.64516 (a/λ), the wavevector is 2.06 (2π µm -1 ). The first-order mode is cut-off, which means that the output silicon waveguide is single mode at this frequency. The dispersion relation ω(k) for index-guided modes with even mirror symmetry in the y-direction was calculated with the use of MPB [START_REF] Johnson | Block-iterative frequency-domain methods for Maxwell's equations in a planewave basis[END_REF]. The gray region represents the light cone (n1= 1.44) and encloses the radiative modes.

Although we treated the problem as a 2D one to reduce computational time and effort, the photonic structure can be extended to a 3D case. A way to treat a 3D structure is to map the 3D problem into an approximative 2D one by using the effective index method [START_REF] Chiang | Dual effective-index method for the analysis of rectangular dielectric waveguides[END_REF]. In this case thus, the index of refraction of the photonic structure is that of the effective index of the fundamental vertical slab mode. For 3D optimizations, we used a 2D optimization domain (i.e., a 2D material grid) and extruded the design into the vertical direction by transferring the optimization domain along all thickness of the silicon layer. The extrusion process consists of stretching the flat, 2D shape vertically to generate the 3D photonic device.

This design method is a convenient method for the optimization and for the fabrication of the device with standard planar fabrication tools such as electron beam lithography and reactive ion etching. We show the results of this approach in the supplementary information. To eliminate the optimization of topological structures with size features unattainable with current nanofabrication tools (i.e., standard electron beam lithography), we include nanofabrication constraints in the design region. Specifically, linear and nonlinear functional transform on the design parameters were included before projecting them onto the simulation domain. A conic filter was used to blur out features with size lower than a threshold with a minimum length of 0.09 µm.

The adjoint solver implemented in Meep is based on the density method and efficiently computes the gradient of an arbitrary function of the electromagnetic fields with respect to ε on a discrete spatial grid.

In the density-based method, each pixel in the design region is either solid or void for high-and low-index materials, respectively. The design is gradually binarized with the application of a series of filters and projectors. One important aspect of this process is that the cost function and its gradient are continuous and thus differentiable. Moreover, the method involves two separate simulations: the first one is the "forward" simulation that computes the fraction of the power emitted by the dipole into the fundamental mode and the fields in the design region and the second simulation is the "adjoint" one that computes the fields in the design region with another special type of source distribution. The gradient of the objective function with respect to the design variables is obtained in a post-processing using the fields in the forward and adjoint simulations. After the calculation of the gradient, we finally run the optimization algorithm based on the method of moving asymptotes (MMA) available in the NLopt library [START_REF] Svanberg | A class of globally convergent optimization methods based on conservative convex separable approximations[END_REF][START_REF] Johnson | The NLopt nonlinear-optimization package[END_REF]. The gradient values are used to update the material density distributions in the design region for the next iteration. Initially, the material distribution within the design region is random and varies within lower and upper bounds from 0 (i.e., 𝜖 1 = 𝑛 1 2 ) to 1 (i.e., 𝜖 2 = 𝑛 2 2 ). The permittivity at each grid cell is a design parameter denoted as 𝜖(𝜌) = 𝜖 1 + 𝜌(𝜖 2 -𝜖 1 ), 𝜌 ∈ [0,1]. Further details on the adjoint method can be found in reference [START_REF] Hammond | High-performance hybrid time/frequency-domain topology optimization for large-scale photonics inverse design[END_REF]. We run all processes for a number of n = 50 runs to statistically treat the modal collection efficiency and LDOS obtained through topology optimization.

We done this by computing the mean m and standard deviation 𝜎 of the objective solutions. The uncertainty is defined as 𝑢 = 𝜎 ⋅ 𝑛 -1 2 (i.e., the standard error).

Results and discussion

To present the results of the optimization algorithm, we divide this section into three subsections, each of one presents the results of a particular source. The comparison between 2D and 3D topology-optimized photonic devices, influence of the design region size and position of the dipole source within the design region are also investigated, specifically in the supplementary materials. Finally, in the supplementary information, we plotted the results for single electric dipole oriented along x and z, and that for a magnetic dipole oriented along the y axis.

Circular polarized dipole source

Figure 2 shows the evolution of the design region as the method finds the best β-factor for one of the topology optimization realizations. The design region follows a rapid transformation from a random initial to a final optimized design. Indeed, we see that the final optimized design reaches its final form from iteration 35 (Fig. 2a). Correspondingly, the β-factor of the topology optimized photonic structure is shown in (Fig. 2b). For collection efficiency we obtain a minimum and maximum of 0.82 and 0.90, respectively. The mean and standard deviation are 0.86 and 0.019. In addition to the best collection efficiency (or β-factor), the optimized photonic structure enhances the LDOS compared to the LDOS in free space at frequency 0.64516(a/λ). For a number of n = 50 realizations, the minimum and maximum LDOS values are 5.8 and 9.7 with a mean of 7.2 and standard deviation of 0.94. The respective uncertainty is 0.13.

In figure 3, we show the distribution of the components of the electric and magnetic field propagating along the optimized structure as well as the absolute squared value of the electric field at the central frequency. We recall here that the algorithm is set to optimize the coupling into the fundamental mode of the output waveguide with out-of-plane magnetic field hy. The circular polarized dipole interacts with the y-even modes, which have non-zero field components ez, hy, and ex. The px and pz components of the dipole couple to the ex and ez components of the available modes. In a regular spin-directional coupling, when the dipole is placed in the near field of a waveguide modes, the interaction is done through the transverse spin angular momentum (SAM) density [START_REF] Bliokh | Extraordinary momentum and spin in evanescent waves[END_REF][START_REF] Bliokh | Transverse and longitudinal angular momenta of light[END_REF][START_REF] Abujetas | Spin angular momentum of guided light induced by transverse confinement and intrinsic helicity[END_REF].

The transverse SAM is either positive or negative inside the waveguide and therefore, by placing the dipole at r0 = (x0, ±z0), the dipole excites the mode propagating along the ±x-axis (Fig. 4b). In optimized design, the topology is formed such that the field modes, specifically the component ex and ez have a fixed amplitude and phase relationship. In the topology-optimized region, the optimization starts first with the dipole in a design region supporting six modes. Because the dipole is placed at the center, the transverse SAM density is zero for all the modes (i.e., the electric part of the SAM Sy = 0). The emission is therefore coupled to both the forward and backward modes. One can see that the first iterations in the optimization generate a mirror that reflects the emission, coupled to the backward modes, to the +x-axis maximizing the collection efficiency.

Because the amplitude and phase relationship between the ex and ez components must be fixed for unidirectional emission, the algorithm enforces this relationship by surrounding the dipole in void and material parts, with some rotational symmetry. The dipole itself has rotational symmetry around the yaxis. The coupling efficiency of a topology optimized structure is three times higher than that with spinorbit locking (see suppl. info.).

To investigate the influence of the size of the design region and the optimization with less available modes in the design region, we conducted optimizations with squared design regions of side 2.0 µm, 1.5µm, and 1.0 µm (i.e., 1.29λ, λ, and 0.64λ). In those cases, the available modes at frequency f = 0.64516 are 4, 3, and 2, respectively. We present in Table II the collection efficiencies into the fundamental mode for the different region side size.

The collection efficiencies decrease as the size of the design region decreases, 12% of the collection efficiency is lost when the side is reduced by a factor of two (i.e., for a region of size of 1.5µm by 1.5µm).

A figure that shows the propagation of the electric field for optimized design regions of 2µm, 1.5µm, and 1µm can be found in suppl. info (Fig. S6). The reduction of the size of the design region to one-wavelength results in a collection efficiency of 73% and a Purcell factor of 10.3. In addition to the influence of the size of the design region, we also investigated the influence of the position of the dipole within the design region on the collection efficiency. Indeed, the topology optimizations in previous cases were done for a dipole placed at the center of the design region. We investigated four different random off-center positions of the dipole in a 3-μm-side design region (Table III). The corresponding electric field distributions and optimized design regions at frequency f = 0.64516 are shown in supplementary information. Although the dipole is placed off-center, the TO leads to high collection efficiency (> 85%) into the fundamental mode while providing values of the Purcell factor comparable to that for the dipole placed at the center of the design region.

For 3D optimizations, we used a 2D optimization domain (i.e., a 2D material grid) and extruded the design into the vertical direction by transferring the optimization domain along all thickness of the silicon layer. The thickness of the silicon layer was fixed to 220 nm and the dipole was placed on top of the design region. The 3D topology optimization with a 2D design region yields collection efficiency as high as 84% and a Purcell factor of 33.6. The collection efficiency is comparable to that obtained in the effective 2D optimization. The Purcell factor increases as the effective mode volume decreases; indeed, the optical confinement is in 3D. The spin-direction coupling via topology optimization is viable in 3D photonic devices.

Huygens' dipole source

Figure 5 shows the topology optimized photonic structure. The minimum and maximum LDOS values are 6.3 and 9.7 with a mean of 7.2 and standard deviation of 0.77. The respective uncertainty is 0.1. For the β-factor we obtained a minimum and a maximum of 0.89 and 0.99, respectively. The mean and standard deviation are 0.98 and 0.04. In figure 5, we show the distribution of the components of the electric and magnetic field propagating along the optimized structure as well as the absolute squared value of the electric field at the central frequency. We recall here that the algorithm is set to optimize the coupling into the fundamental mode of the output waveguide with out-of-plane magnetic field hy. The Huygens dipole produces a net power flow in each direction due to the relative amplitude and phase between the electric and magnetic dipoles. Here, the electric and magnetic dipoles are orthogonal and inphase, which produces directionality when placed in the near field of a waveguide (Fig. S12b andd).

The coupling efficiency of a Huygens dipole placed at the center of the topology optimized region is more than one order of magnitude higher than that of the Huygens dipole placed in the near field of a waveguide, with an enhancement of about 15 (Fig. S12f).

In the Huygens dipole, the pz and hy components of the dipole couple to the transversal ez and hy components of the available modes. It does not involve coupling to the longitudinal x-component of the electric field. In optimized design, the topology is formed such that the field modes, specifically the component ez and hy have a fixed amplitude and phase relationship. The emission is therefore coupled to the forward modes but some of the radiation radiates at some angle. One can think that optimization generates a focusing structure to focus the light into the output waveguide, maximizing the collection efficiency. Because the amplitude and phase relationship between the ez and hy components must be fixed, the algorithm enforces this relationship by surrounding the dipole in material part, with some mirror-symmetry with respect to the z = 0 plane, as the Huygens dipole itself.

Finally, we conducted the effective 2D optimization (2D TO then 3D simulation with the geometry extruded in the vertical direction) with the Huygens dipole on the top of the optimization domain. We obtained a coupling efficiency of 72% and a Purcell factor of 22.8.

Janus dipole source: coupling and noncoupling faces

In table IV, we present the basic statistics for n = 50 realization of the Janus dipole source. Figure 7 shows the topology optimized photonic structure. We recall here that the algorithm is set to optimize the coupling into the fundamental mode of the output waveguide with out-of-plane magnetic field hy. The Janus dipole has orthogonal electric and magnetic dipoles with a π/2 rad phase shift. According to the positive or negative phase shift, the emission couples either into the radiative modes (i.e., Janus noncoupling) or into the guided modes (i.e., Janus coupling) when placed near a waveguide. We calculated that the coupling efficiency of a coupling Janus dipole placed at the center of the topology optimized region is more than one order of magnitude higher than that of the coupling Janus dipole placed in the near field of a waveguide, with an enhancement of about 25 at frequency 0.64516(a/λ) (Fig. S15h). For the noncoupling Janus dipole, the enhancement is more than 1200 times. The topology optimization strategy produces the region that allows a noncoupling Janus dipole to couple its emission into the guided mode. Similarly, in the Janus dipole, the px and hy components of the dipole couple to the ex and hy components of the available modes. In optimized design, the topology is formed such that the field modes, specifically the component ex and hy have a fixed amplitude and phase relationship. The emission is therefore coupled to the forward and backward modes. Also, the optimization generates a focusing structure on the output waveguide, maximizing the collection efficiency. A collection efficiency of 70% was obtained for a 3D simulation with the effective 2D optimization (2D TO then 3D simulation with the geometry extruded in the vertical direction) and a Purcell factor of 17.

Conclusion

Topology optimization strategies were used to find the best geometries that can increase not only the coupling efficiency of dipole emission radiation into a desired photonic mode but also to enhance the local density of states (LDOS). Linear and circular polarized dipoles placed in small topology-optimized photonic structures (i.e., design region of 3 µm x 3 µm) were used. We demonstrated a mean coupling efficiency and LDOS enhancement higher than 0.86 and 7.2, respectively, for circular polarized, Huygens, and Janus dipoles. In addition, the coupling efficiency of circular polarized, Huygens, and Janus dipoles placed at the center of a topology optimized region is more than one order of magnitude higher than that of the same dipoles placed in the near field of a waveguide. The approach adopted to design 3D photonic device by the extrusion of the optimization design into the vertical direction simplifies the optimization process for fabrication with planar technologies. For practical applications, it is important to use quantum emitters with physical size lower than the volume of a nearby optical field mode in the optimized photonic structure to consider the emitter as a point-like dipole. We expect to extend our approach to other integrated photonic material platforms such as silicon nitride and silicon carbide integrated photonics.

Figure 1 .

 1 Figure 1. Topology optimization problem. (a) 3D schematic of the photonic structure with a square optimization region. (b) The 2D problem geometry where the 3 µm by 3 µm optimization region is seen. The right side of the optimization region is connected to the output port waveguide where the power coupled into the fundamental mode is calculated. The red profile is the calculated eigenmode at frequency 0.64516 (a/λ). The wavevector at this frequency is 2.06 (2π µm -1 ). The dipole is represented by a red point and placed at the center of the optimization region. The light green box represents the flux monitor to calculate the power radiated by the dipole. The green region around the structure represents the PML. (c) Calculated dispersion relation of the output silicon waveguide.

Figure 2 .

 2 Figure 2. Topology optimization of the integrated circular polarized single photon source. (a) Evolution of the design region where the position r0 of the dipole is at the center of the 3 µm by 3 µm design region. (b) Scatter plot with circular markers at the locations specified by the enhancement of the local density of states (i.e., the Purcell factor) and the collection efficiency at frequency 0.64516(a/λ) for 50 optimization runs. Each optimization is done for a broad band spectrum centered at frequency 0.64516(a/λ) and with a width of 0.129. Inset: the collection efficiency and Purcell factor spectra for the optimization run number of 20.

Figure 3 .

 3 Figure 3. Distribution of the field along the topology optimized integrated single photon source. Propagation of (a) the x-component, (b) z-component of the electric field and (c) the y-component of the magnetic field along the optimized structure. (d) Distribution of the electric field squared. The scale bar is 500 nm in (a) and the inset in (b) shows the orientation of the dipoles. The frequency is 0.64516(a/λ).

Figure 4 .

 4 Figure 4. Eigenmode analysis and FDTD calculation of the spin-orbit coupling. (a) Distribution of the ratio between ie x and e z field components. The points in ie x /e z show the region where the amplitude of these two components is equal at z0 = ±190 nm. (b) Calculated distribution of the only nonzero components of the spin AM density Se (i.e., Sy, the y-component, green solid line) and the spin momentum Psx. Distribution of the electric field for a circular polarized dipole placed at (c) z0 = 190 nm and (d) z0 = 0 in a straight waveguide, and (e) placed at the center of the topology optimized. Unidirectional coupling into the forward mode is obtained at the spin-orbit locking position in (c) and enhanced unidirectionality in the optimized structure (i.e., in (e)). The distributions are obtained at frequency 0.64516(a/λ) and the scale bar is 1 µm.

Figure 5 .

 5 Figure 5. Topology optimization of the integrated Huygens' dipole. (a) Evolution of the design region. (b) Enhancement of the local density of states and collection efficiency corresponding to the optimized photonic structure on one of the topology optimization realizations. Inset: Purcell and collection efficiency spectra.

Figure 6 .

 6 Figure 6. Distribution of the fields along the topology optimized integrated Huygens' dipole. Propagation of (a) the x-and (b) z-component of the electric field, real part. (c) Distribution of the y-component of the magnetic field and (d) squared electric field. The scale bar is 500 nm and the inset in (b) shows the orientation of the dipole. The frequency is 0.64516(a/λ).

Figure 7 .

 7 Figure 7. Topology optimization of the integrated Janus single photon source. Evolution of the design region for (a) a noncoupling and (b) coupling Janus dipole. (c) Enhancement of the local density of states and collection efficiency corresponding to the optimized photonic structure on one of the topology optimization realizations.

Figure 8 .

 8 Figure 8. Topology optimization of an integrated Janus single photon source. Purcell-and β-factors spectra at run n = 20 for Janus dipole (a) noncoupling and (b) coupling. Distribution of the electric fields along the topology optimized integrated source. Propagation of the real part of the y-component of the Janus dipole (c) noncoupling and (d) coupling. The inset in (c) and (d) shows the orientation of the Janus dipoles. The frequency is 0.64516(a/λ).

Table II .

 II Collection efficiencies and Purcell factors for topology optimized regions with different sides.

	Design region side (µm)	Collection efficiency	Purcell factor
	3.0	0.88	9.4
	2.0	0.77	7.4
	1.5	0.76	5.9
	1.0	0.73	10.3

Table III .

 III Collection efficiency at f=0.64516 in units of (a/λ) for four different positions around the design region of the circular polarized dipole. Δx and Δz are the offset distance relative to the center of the design region.

	Offset from the center	Δx (µm)	Δz (µm)	Collection efficiency	Purcell factor
	0	0	0	0.88	10.1
	1	-0.75	+0.75	0.85	8.2
	2	+0.75	+0.75	0.92	9.2
	3	-0.5	-1.0	0.94	7.8
	4	+1.25	-0.5	0.90	2.65

Table IV .

 IV Basic statistics on n = 50 optimization runs for a Janus dipole source.

			noncoupling		coupling
		LDOS Collection	LDOS	Collection
			efficiency		efficiency
	min	5.9	0.91	6.9	0.85
	max	15.6	1.0	14.6	1.0
	mean	8.1	0.98	10.0	0.92
	Std. deviation	1.7	0.03	1.85	0.04
	uncertainty	0.24	0.005	0.26	0.005

It is interesting to note that the optimization algorithm works for the noncoupling face Janus dipole. The noncoupling face presents values comparable to that of the coupling Janus dipole.

Acknowledgments

We conducted the calculation with the resources of the French regional Grand-Est HPC Center ROMEO.

Funding

We thank the China Scholarship Council File No. 202 008 070 054 for partial financial support. We conducted this work within the frame of the Graduate School (Ecole Universitaire de Recherche) "NANO-PHOT," contract ANR-18-EURE-0013.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Disclosures

The authors declare no conflicts of interest.

Supplemental document

Dispersion relation in a 3-µm-thick vertical slab waveguide. Intrinsic properties of eigenmodes in waveguides. Comparison between non-optimized and topology optimized regions. Influence of the design region size and dipole position on the collection efficiency for the circularly polarized dipole. Threedimensional (3D) topology-optimized photonic device. More results for individual electric and magnetic dipoles oriented along three different axis l = x, y, and z.