Numerical optimization of the extraction efficiency of a quantum-dot based single-photon emitter into a single-mode fiber

We present a numerical method for the accurate and efficient simulation of strongly localized light sources, such as quantum dots, embedded in dielectric micro-optical structures. We apply the method in order to optimize the photon extraction efficiency of a single-photon emitter consisting of a quantum dot embedded into a multi-layer stack with further lateral structures. Furthermore, we present methods to study the robustness of the extraction efficiency with respect to fabrication errors and defects.


Introduction
Single-photon emitters (SPE) are essential building blocks of future photonic and quantum optical devices. Solid-state SPEs, such as self-assembled quantum dots (QDs) and defect centers in solids, provide a scalable platform with outstanding optical properties in terms of the suppression of multiphoton emission and photon indistinguishability [1]. However, not only the properties of the SPEs but also the solid-state structures surrounding the QDs have an important influence on the system's performance. The surrounding structure can enhance or suppress the rate of spontaneous emission, a process known as Purcell effect [2]. Moreover, the efficiencies with which emitted photons are extracted into a specific direction or coupled into an optical fiber depend strongly on the geometry of the surrounding structure [3].
Advanced nanoprocessing technologies allow for the realization of geometries with many degrees of freedom. However, from an experimental and technological perspective it is often not clear which geometries and geometric parameters lead to optimal results. In this case the numerical simulation of the micro-optical structures can give important insights [15][16][17][18].
In the following we present a method to simulate the light field emitted by a QD embedded into micro-optical structures and its coupling efficiency to an external optical fiber. The coupling efficiency is very sensitive to the specific intensity and phase profile of the emitted light field when it enters the optical fiber. Hence, the numerical simulations require a high accuracy [19]. The presented numerical approach is based on the finite element method (FEM) in the frequency domain [20]. For nano-optical resonant devices, this method offers several important benefits compared to alternative methods, like Finite Difference Time Domain (FDTD) or rigorous coupled-wave analysis (RCWA), with respect to accuracy and computation time. Firstly, the geometry is discretized with a non-uniform mesh allowing for an efficient and highly accurate modeling of the material interfaces without stair-casing effects as for FDTD and RCWA [21]. Secondly, the field is represented with local polynomial ansatz functions of various polynomial degree on the patches of the mesh. In this way and in contrast to standard FDTD, highly non-uniform field profiles can be captured by locally refined meshes, whereas regions of wave propagation can be very efficiently discretized using high-degree polynomial ansatz functions (p-refinement) [22]. Thirdly, since we are only interested in a QD emitting at a specific frequency, it is advantageous to solve Maxwell's equation directly in frequency domain. This is again in contrast to FDTD where the frequency response is reconstructed by a spectral analysis of a broadband excitation, which leads to long computation times especially for fine meshes required for high accuracies. As a drawback, FEM has a higher memory footprint compared to FDTD. However, this was not the limiting factor for the studied SPE. In order to further improve the numerical convergence of the FEM approach, we exploit the radial symmetry of the system and use a subtraction method to cope with the singular behaviour of dipolar field emitted by the QD.
We apply the numerical method in order to study QDs embedded into micro-optical structures, which are fabricated on top of a Bragg reflector. The objective of the simulations is to identify geometries that efficiently couple the emitted light of the QD into an optical fiber above the structure to further improve performance of stand-alone fiber-coupled single photon sources [23] in the future. For clarity reasons, we focus on two embedding structure designs, microlenses and mesas, which can be produced with high accuracy utilizing 3D in-situ electron-beam lithography [15,24].
The paper is organized as follows. After the introduction of the physical system of the embedded QD in Sec. 2 we present the numerical approach in Sec. 3. The approach is based on a finite element method (FEM) [25] and takes the highly singular nature of the QD light-field into account. We demonstrate that the computational times can be significantly reduced by exploiting the cylindrical symmetry of the considered geometry. In Sec. 4 we apply the method in order to optimize the system parameters. By scanning the parameter space around the optima we assess the robustness of the geometry against the influence of fabrication errors.

System description
In the following we consider a system consisting of a single QD emitting at a vacuum wavelength of λ = 1, 300 nm in the telecom O-band. The QD is embedded into a spherical lens or a mesa structure made from gallium arsenide (GaAs; refractive index n struct = 3.4). An underlying Bragg reflector made from layers of GaAs (n Bragg1 = 3.4) and aluminum gallium arsenide (Al 0.9 Ga 0.1 As; n Bragg2 = 3.0) reflects the light emitted by the QD back into the upper hemisphere (n space = 1.0). The light is coupled into an optical fiber above the QD consisting of a homogeneous fiber core and a homogeneous fiber cladding (n core = 1.5, n clad = 1.45, NA ≡ n 2 core − n 2 clad = 0.38). Figure 1 shows the two different types of diffractive structures -spherical lens and mesa. The aim of these structures is to direct as much light as possible into the optical fiber. In Sec. 4 the system behavior will be studied as a function of five geometry parameters (see Fig. 1). These include the fiber-core diameter and the distance between fiber and spherical lens or mesa as well as the geometrical parameters of the embedding structure (dimensions of the spherical lens or mesa and the elevation of the QD). The geometry of the embedding structure can be controlled through the growth process and 3D in-situ electron-beam lithography [15].
In the following, we consider fiber-core diameters between 1, 000 nm and 2, 500 nm. Optical fibers with such small core diameters can be fabricated, e. g., by flame heating and pulling standard fibers into a narrow thread [26]. The mode field diameter (MFD) of the single-mode fiber, i. e. the diameter where the field energy density has dropped to 1/e of the maximal field energy density, depends on the fiber-core diameter. For the considered parameter range the MFD varies between 2, 400 nm and 3, 400 nm. The simulation results in Sec. 4 show that the relatively small fiber-core diameters and correspondingly small MFDs are favourable in order to gain a large overlap between the light field emitted by the QD and the propagation mode of the fiber. The considered system consists of a QD dipole source (red point) embedded into a diffractive structure (GaAs, blue), a Bragg reflector (alternating layers of GaAs (blue) and Al 0.9 Ga 0.1 As (gray)), and an optical fiber with homogeneous fiber core (orange) and fiber cladding (yellow). The Bragg reflector is grown on a substrate made of GaAs (blue) and has a GaAs top layer (blue, thickness = 195 nm ). The system is parametrized by 5 length scales, the diameter of the fiber core (d core ), the width and height of the lens or mesa (w lens , h lens ), the elevation of the dipole (h dip ), and the distance between lens and fiber (s lf ). Left: Spherical lens setup. Right: Mesa setup.

Approximation of the quantum dot by a classical dipole
The self-assembled QD under consideration consists of InGaAs and has a typical extension of 20 nm in horizontal direction and 5 nm in vertical direction. After excitation, it traps a single electron-hole pair, which can recombine through the spontaneous emission of a photon. To treat the interaction between the QD and the light field in the diffractive structure rigorously, one would need to invoke quantum electrodynamics (QED) [27]. However, the structure into which the photon is emitted has a low quality factor. That is, the photon energy diffuses quickly without having the time to re-interact with the QD. In this case one can model the emission properties of the QD by an oscillating point-like current density Here, r QD is the position of the QD and ω its emission frequency. The QD has a larger lateral than vertical extension. Hence, the electronic state is excited horizontally and the dipole moment p lies in the horizontal plane. A more detailed explanation of the dipole approximation is given in Appendix A.

Time-harmonic Maxwell's equations
Since the considered source current density J(r, t) is time-harmonic (see Eq. (1)), the same holds for the electromagnetic fields E(r, t) and H(r, t). Therefore, we consider the Maxwell's equations for the electric field in the frequency domain [21] The permeability µ = µ(r, ω) and the permittivity = (r, ω) depend on the material properties and the considered wavelength λ = 2π/ω of the light emitted by the QD. For simplicity we will henceforth implicitly assume that the fields and material tensors depend on ω. The goal of the numerical method is to find the field distribution E(r) for different geometrical parameters of the system, and to determine in a second step the coupling efficiency of the field into the optical fiber.

Subtraction method
The electric field produced by the dipole emitter defined in Eq. (1) diverges for r → r QD . Hence, a direct finite-element discretization of E suffers from a slow numerical convergence. To cure this, the electric field is expressed as a sum E = E s + E c of the singular dipole field E s and a correction field E c [28]. The singular field is chosen to be the analytically known dipole-field solution of the Maxwell's equations with constant material tensors µ d = µ(r QD ) and d = (r QD ) . As shown in [28], the Maxwell's equation for the correction field reads In the vicinity of the dipole emitter, the right-hand side of Eq. (4) vanishes since µ → µ d and → d . Therefore, the correction field is smooth close to the dipole and can be efficiently computed using the finite-element method decribed in [25].

Symmetry adaption
For the considered setup, the material tensors are rotationally symmetric, i. e. in cylindrical coordinates it holds µ(r, z, φ) = µ(r, z) and (r, z, φ) = (r, z) with r = |r|. Therefore, we expand the electric field in a series of Fourier modes Likewise, the right hand side of Eq. (4), can be expanded into a series Inserting Eqs. (5) and (6) into Eq. (4) and integrating over 1 2π ∫ dφ e −imφ yields independent 2-dimensional differential equations for every Fourier mode m ∈ Z, i. e. 1 2π where ∇ arises from ∇ by the replacement ∂ φ → in. Due to the reduced dimensionality of the equations the corresponding finite-element computations require much less computation time and reach a higher accuracy level than in the case Maxwell's equations are solved on a 3D mesh (see Sec. 3.6).
In numerical computations, only a finite number of Fourier modes is computed. The number of Fourier modes is chosen automatically by an adaptive algorithm such that the fields converge within some level of accuracy. For the case that the dipole is positioned at the symmetry axis (r QD = (0, 0, z) T ), only the Fourier modes m = −1, 0, 1 contribute to the electric field expansion. However, we like to stress that even though the expansion relies on a cylindrical symmetric geometry, the source currents and light fields do not need to be symmetric. For example, also off-axis dipoles can be simulated by the method (see section 4).

Coupling efficiency to fiber modes
The geometry of the optical fiber as described in section 2 is invariant in z-direction. The eigenmodes E n (r) for a given frequency ω are solutions of Eq. (2) with no source current, i. e.
For practical reasons one is interested in the fraction of the light field energy that is coupled into the fundamental mode of the optical fiber. In order to define the coupling efficiency, we expand the electromagnetic field, which is scattered into the fiber, into a series of eigenmodes, i. e., As scalar product we employ the following overlap integral: where the integration is performed over a cross section of the optical fiber.
In Appendix B we show that the modes are orthonormal with respect to the scalar product of Eq. (11). Due to the orthonormality, the total power P = 1 2 Re ∫ n · (E scatt × H * scatt ) emitted into the optical fiber is given as In the case of guided modes in a loss-free non-active medium the eigenmodes E n are real-valued and the expression simplifies to The coupling efficiency towards a specific mode n is defined as the powerflux P n = Re E n , E scatt E n , E * scatt to this mode divided by the total emitted power of the dipole P tot . In the following the coupling efficiency to the two degenerate fundamental eigenmodes E 1 and E 2 of the optical fiber is considered. We would like to note, that the equations are not only valid for optical fibers, but for any waveguide structure.

Validation and convergence
In order to validate the presented methods, we compare results of the symmetry-adapted simulations in effectively two dimensions (2D) with those of non-symmetry-adapted threedimensional simulations (3D). To this end we consider a simplified setup consisting of a Bragg mirror with 10 double-layers (n Bragg1 = 3.4, n Bragg2 = 3.0). The diffractive structure (spherical lens or mesa) is replaced by a flat GaAs layer with a thickness of 395 nm (n struct = 3.4). The dipole is placed at at a height of 195 nm within this layer and the fiber with core diameter 1, 500 nm (n core = 1.5, n clad = 1.45) has a distance of 100 nm to the flat GaAs layer (n space = 1.5). This setup has a coupling efficiency of η = 3.8% to the two degenerate fundamental modes of the fiber. Figure 2 shows the convergence behavior for an increase of the polynomial degree p of the finite-element ansatz functions (finite-element degree) [22] for the 2D and 3D computations. Both computations converge to high numerical accuracies on the order of an absolute error of 10 −4 showing that the methods provide reliable numerical results for the coupling efficiency. In order to determine the coupling efficiency with an accuracy better than 1% a polynomial degree of p = 3 (2D case) or p = 4 (3D case) is sufficient. For this accuracy level the computation times of the 2D and 3D case differ drastically. While the 3D calculations require about 1 hour with 8 CPU threads, the 2D computation requires less than 10 seconds. We attribute the slower convergence of the 3D calculations to the fact that the meshing of the structure in the 3D calculations leads to an additional discretisation error along the angular coordinate φ.
In the next section, the speed-up due to the exploitation of the cylindrical symmetry will enable the exploration of the behavior of the coupling efficiency for a large number of different parameter values of the geometry.

Numerical optimization of the geometry
In the following, we seek to optimize the geometry of the micro-optical structures such that the light emitted by the QD is efficiently coupled into the two fundamental modes of the fiber. The free parameters under consideration are the width w lens and height h lens of the diffractive structures (spherical lens or mesa), the elevation of the dipole h dip within these structures, the distance between lens and fiber s lf , and the diameter of the fiber core d core (see figure 1). The numerical computations of the spherical lens and mesa geometries where performed on a computational domain with a radius of 6, 000 nm and transparent boundary conditions. The systems were discretized with about 3, 700 triangular mesh elements and required computation times of about 40 seconds.
The short computation times of the 2D calculations allow for a rather fine sampling of the parameter space. However, a sampling on a 5-dimensional grid would still require too many computations. Therefore, we first sample the parameters of the lower part of the setup (w lens , h lens , and h dip ) for fixed parameters d core = 1, 500 nm and s lf = 100 nm. The results of the parameter scans are shown in Fig. 3. Obviously, the coupling efficiency depends very sensitively on the width and height of the spherical lens and mesa. Especially the complex behavior for a variation of the width of the mesa is an indication that several resonant states within the structure can be exited, each having a different overlap with the fundamental propagation modes of the optical fiber. In this step the best coupling efficiencies obtained are η = 26.2% for the spherical lens setup and η = 21.0% for the mesa setup.
In the second step, the diameter of the fiber core d core and the distance between lens and fiber s lf are sampled while setting the values of w lens , h lens , and h dip to the values optimized in the first step. The results of the parameter scans are shown in Fig. 4. The spherical lens setup is more sensitive on the lens-fiber distance than the mesa setup, indicating that the field radiated upwards has a stronger divergent behavior. For large fiber-core diameters the coupling efficiency of both structures slowly decreases showing that narrow tapered fibers are crucial in order to reach high coupling efficiencies. After the second optimization step the best coupling efficiencies obtained are η = 29.9% for the spherical lens setup and η = 23.2% for the mesa setup. The optimized parameters are summarized in Tab. 1. In Fig. 5 the energy density distribution is visualized for the optimal geometries of the spherical lens and the mesa setups. For the optimal mesa setup less radiation is propagating into the Bragg reflector. Nevertheless, the spherical lens leads to a considerably larger coupling efficiency of about 30% as compared to 23% for the setup with the mesa. This indicates that the specific intensity and phase profile of the light field entering the fiber is of importance.
Although the spherical lens setup is performing better, the structure is much harder to produce  Fig. 5. Visualization of the energy density on a logarithmic scale of the light field emitted by the QD embedded in the optimized spherical lens (left) and mesa (right). A cut through the optimized geometry is shown in front of the energy-density plots with the fiber core in purple, the fiber cladding in red, the spherical lens or mesa structure in green and the layers of the Bragg reflector in light and dark blue. The optimal mesa has a smaller coupling efficiency although less radiation is propagating into the Bragg reflector. This indicates that the specific intensity and phase profile of the light field entering the fiber is of importance. than a mesa. Moreover, small errors in the fabrication process can have a large influence on the geometry of the spherical lens. Therefore, in the next section we study the sensitivity of the coupling efficiency with respect to fabrication errors.

Sensitivity analysis
The fabrication of the diffractive structures is subject to various errors. While the heights of the epitactically grown structures can be produced with a high accuracy of a few nanometers, the lateral accuracy is generally lower. For example, the lateral position of the dipole has an uncertainty of about 34 nm [29] and the widths of the spherical lens and mesa can be controlled with an accuracy of about 100 nm. Moreover, also the lens-fiber distance has a relatively large uncertainty of about 50 nm. The sensitivities with respect to errors in the system parameters are depicted in Fig. 6. The figure shows that the spherical-lens setup is rather robust against fabrication errors of the lens width, the lens-fiber distance, the fiber-core diameter, and the lateral mismatch of the QD. On the other hand, an error in the lens height or the dipole elevation leads to a severe degradation of the coupling efficiency.
In comparison, the coupling efficiency of the optimized mesa structure is more strongly influenced by fabrication errors of the width of the mesa and the lateral position of the QD (see right-hand side of Fig. 6). Also in this case, errors in the lens-fiber distance and the fiber-core diameter hardly influence the coupling efficiency. Altogether, both structures are relatively robust against lateral fabrication errors and thus meet current processing accuracy tolerances. Noteworthy, the right-hand side of Fig. 6 reveals, that the simple two-step optimization scheme did not converge to a local maximum. After the second optimization step (variation of the fiber-core diameter and the lens-fiber distance), the optimized values of the first step are not at the local maximum any more. For example, a lower value of the mesa width leads to a small increase of the coupling efficiency. Global optimization schemes, such as Bayesian optimization, tackle these problems by varying all system parameters over a given domain [30] and can thus further enhance the achievable extraction efficiency. However, the application of such methods is beyond the scope of this study.
The fabricated structures do not only suffer from errors of the geometry parameters, but also from an etching-induced surface roughness. Especially, the fabrication of spherical lenses can lead to non-negligible defects of the lens surface. The surface roughness breaks the cylindrical symmetry such that the corresponding effects can only be assessed by performing full 3D simulations. Therefore, we consider a reduced geometry, consisting of a Bragg reflector with only 10 instead of 25 double-layers. Although this reduces the optimal coupling efficiency to 19%, the relative magnitude of roughness effects should be comparable between the full setup and the reduced one.
We assume that the surface roughness has a Gaussian autocorrelation. Hence, it can be parametrized by a single correlation length and a roughness amplitude [31]. For the numerical study we vary the roughness amplitude between 0 and 50 nm and set the correlation length to l = 200 nm. The mean influence and the variance are estimated by generating 10 random roughness profiles for each amplitude. Figure 7 shows how the coupling efficiency degrades with increasing roughness amplitude. With increasing amplitude also the variance of the degradation increases. The numerical results show that surface roughnesses with an amplitude below 10 nm can be usually tolerated, while larger amplitudes quickly lead to a severe drop of the coupling efficiency.

Conclusion
We have presented an efficient numerical method for the simulation of localized light sources embedded in micro-optical structures, taking into account both, the highly singular nature of the emitted light field as well as the rotational symmetry of the geometrical setup. The method was applied for studying the behavior of self-assembled quantum dots embedded into two types of diffractive structures -spherical lens and mesa -placed above a Bragg reflector. We have performed a sampling of the parameter space of the system geometry and an optimization of the structure with respect to the coupling efficiency to the fundamental modes of a fiber. Moreover, we studied the robustness of the optimized geometries in the presence of fabrication errors. The optimized geometries for the spherical lens and mesa setups presented here are not necessarily the global optima of the parameter space. Although the computation times are relatively small, a high-resolved grid search of the optimum in the full parameter space is still unfeasible. Global optimization schemes, such as Bayesian optimization, tackle this problem by sampling the parameter space efficiently based on a statistical model of the objective function [30]. For the future we plan to apply Bayesian optimization methods in order to extend the covered parameter space and to identify geometries with an even higher coupling efficiency.

Appendix A -dipole approximation
In Coulomb gauge, the interaction between a charge distribution of a QD and the light field is determined by the HamiltonianĤ where A is the vector potential of the electromagnetic field and q n , m n and r n are the charge, mass and position of the charges, respectively. The dynamics of the system is governed by the coupling matrix elements of this operator in the basis of the Fock states of the electromagnetic field |n and the ground and excited state Ψ g and |Ψ e of the QD. Within the dipole approximation one neglects the variation of A(r) over the extend of the charge distribution. That is, in the integration over the electronic degrees of freedom one replaces A(r) by A(r QD ), where r QD is the position of the QD. In this case the properties of the QD are described by its emission frequency ω and its dipole moment Since the QD considered in this work has a larger horizontal than lateral extension the first state |Ψ e is excited in a horizontal direction and the dipole moment in z-direction can be neglected. The dipole approximation is only valid if gradients of the electric field can be neglected around the QD. Therefore, the wavelength in the embedding material λ/n struct and the length scales of the surrounding geometry must be large compared to the dimensions of the QD. For a QD that is emitting with a vacuum wavelength of 1, 300 nm we have λ/n struct = 382 nm, which is an order of magnitude larger than the dimensions of the QDs (30 nm). Moreover, the comparison with experimental data shows that the dipole approximation provides reliable results as long as a QD is separated by more than roughly 100 nm from a material interface [32,33].
If the QD is placed inside an optical cavity, the coupling strength between the QD and the light field can be described by the coupling constant where V is the volume of the cavity [27]. The strong coupling regime satisfies the condition α γ, where γ is the photon decay rate inside the cavity. In this case the dynamical behavior of the quantum dot is strongly influenced by the surrounding cavity, such that only a QED treatment can accurately describe the system dynamics. A strong coupling leads, for example, to a frequency splitting of the emission spectrum. This effect cannot be reproduced by replacing the QD with a classical light source.
For the optical structures considered in this study, the photon energy quickly radiates into the upper hemisphere. For this weak coupling regime (α γ) it can be shown that QED and the classical theory give the same results for the behavior of the spontaneous emission in the cavity [27]. In this case the dipole can be described by a classical oscillating dipole moment µ(t) = Re pe −iωt e −γt/2 .
The decay rate γ inside a cavity is typically in the GHz regime. At optical frequencies the effects of the decay can be neglected, such that the QD can be modeled by an oscillating dipole with the point-like current density

Appendix B -waveguide mode orthonormality
In the following we will show, that the eigenmodes of a waveguide structure, such as an optical fiber, are orthonormal with respect to the scalar product where the integration is performed over the cross section of the waveguide (w. l. o. g. at z = 0). The eigenmodes E n (r) for a given frequency ω are solutions of Eq. (2) with no source current, i. e. ∇ × µ −1 ∇ × E n (r) − ω 2 E n (r) = 0.
It can be easily confirmed that for any solution E n also the field reflected at the x − y-plane at z = 0 is a solution of of Eq. (20).
We consider the integral of the product of a reflected eigenmode E m and the left-hand side of Eq. (20) over the volume along the waveguide from z = 0 to z = δ