Ab initio electronic stopping power for protons in Ga0.5In0.5P/GaAs/Ge triple-junction solar cells for space applications

Motivated by the radiation damage of solar panels in space, firstly, the results of Monte Carlo particle transport simulations are presented for proton impact on triple-junction Ga0.5In0.5P/GaAs/Ge solar cells, showing the proton projectile penetration in the cells as a function of energy. It is followed by a systematic ab initio investigation of the electronic stopping power (ESP) for protons in different layers of the cell at the relevant velocities via real-time time-dependent density functional theory calculations. The ESP is found to depend significantly on different channelling conditions, which should affect the low-velocity damage predictions, and which are understood in terms of impact parameter and electron density along the path. Additionally, we explore the effect of the interface between the layers of the multilayer structure on the energy loss of a proton, along with the effect of strain in the lattice-matched solar cell. Both effects are found to be small compared with the main bulk effect. The interface energy loss has been found to increase with decreasing proton velocity, and in one case, there is an effective interface energy gain.


Introduction
With the new advancements in solar cell technology and the rapid expansion of space missions, understanding the key aspects which link the macroscopic response of solar cells of current and future spacecraft to the fundamental processes of particle stopping inside the target has acquired new relevance.Solar cells are key components of spacecraft and satellites as they provide power for navigation, communication, data handling, thermal control, and functioning of the instrumentation.At present, the so-called triple-junction (TJ) solar cells, made of Ga 0.5 In 0.5 P, GaAs, and Ge layers (see scheme in figure 1), are the state of the art for space applications, generally protected by an amorphous SiO 2 coverglass.All layers are usually lattice-matched to the substrate (Ge).
In space, radiation degrades the performance of the solar cells mainly by cumulative effects involving atomic displacements.Such displacements introduce consequent defect levels in the electronic structure, which in turn affect the output current.The atomic displacements are caused by the scattering of both primary incoming charges (protons) and secondary projectiles (displaced target atoms, protons, neutrons) generated by the primary radiation.

µm GaAs
Figure 1: Schematic structure of the triple-junction solar cell.The active layers (GaAs and Ga 0.5 In 0.5 P) are lattice matched to the Ge substrate.The tunneling junctions are usually made of p-and n-doped GaAs.Metallic contacts are often aluminum.A detailed structure can be found in Ref [1].
energies below 10-30 MeV (depending on the material) the damaging interaction is mainly due to the Coulomb repulsion between the proton and the target without affecting the nuclei, while hadronic interactions, involving sub-nuclear particles, take over at energies above 30-50 MeV.The Coulomb contribution to the non-ionizing energy loss (NIEL, the rate of energy deposition to atomic displacements, linked to the nuclear stopping) is the main one in the relevant ground testing energy range (0-10 MeV), using unidirectional fluxes on unshielded solar cells.Such energy range and configuration have been shown to be representative of the cumulative damage caused by omnidirectional space radiation on shielded solar cells, where energies in general go up to hundreds of MeV [2,3,4].
The (Coulomb) NIEL is often extracted from SRIM (The Stopping and Range of Ions in Matter) [5,6], from tabulateddata-based approaches like SR-NIEL [7], or from Monte Carlo particle transport codes such as Geant4 [8,9,10] (the latter two also provide the hadronic contribution).Despite some differences, some basic assumptions are similar in these approaches: crystalline order is not considered for the target (thus, there is no distinction between channeling and off-channeling conditions), cascades are treated within the binary-collision approximation (BCA), and the description of electronic stopping is based on the (perturbative) Bethe-Bloch theory at high energy and by the Lindhard dielectric theory combined with the local-density approximation (LDA) at lower energies (E/n < 2 MeV), implying that the electronic stopping power (ESP) at each point in the system is considered to be the same as that of a uniform electron gas of the same density.An ad-hoc parametrization of Lindhard results is used at intermediate energies.In general, the low energy physics in both SRIM and Monte Carlo particle transport codes, such as Geant4, strongly relies on phenomenology and the use of experimental data.
The influence of the ESP on defect production by irradiation of different ions at different energies is under intense investigation [11,12,13,14,15,16,17,18,19,20,21,22,23,24].Recent studies, which address the ESP problem from the ab initio point of view and combine it with the molecular dynamics simulations of damage cascades, show the importance of including electronic stopping effects in cascade formation, in terms of the number of formed defects and cascade morphology [22,23,24].Electronic excitation, as a result of the moving primary and secondary projectiles in the first stages of the radiation damage process, can significantly alter the subsequent ionic displacements and the number of defects formed in the solid.Thus, accurate description of the ESP processes is essential for the correct estimation of the radiation effect in materials.
Real-time time-dependent density-functional theory (RT-TDDFT) has opened the way for the predictive first-principles description of electronic stopping processes in condensed-matter systems, in a non-linear and non-perturbative way [25,27,28,29,30,31,32].Recent RT-TDDFT studies on channeling conditions for proton impact have found a strong sensitivity of the ESP on crystal direction and impact parameter [33].Deviations from velocity proportionality at very low velocity have also been obtained [33,34], consistent with the finite-velocity threshold observed in experiments [35].Others have highlighted the different contributions of semicore electrons in channeling versus offchanneling conditions [36].Local enhancement and reduction of the ESP for specific channels and velocity-dependent deviations from an ideal channeling trajectory were also reported [32], as well as the influence of a strongly anisotropic crystal structure on the stopping power [37].
The explicit consideration of channeling trajectories in the computer simulation of the radiation damage has been shown to be important [38].Both experiments and simulations show that ion channeling plays a significant role in radiation damage cascades in crystalline materials.It has been observed that for an ion moving with energies in the keV range, the penetration depth is much larger due to the channeling effect, which affects the spatial distribution and the shape of the displacement cascade [39,40,41].The channeling effect also contributes to the formation of the sub-cascades [42,43,41] and increases the time of the defect formation [44].
An initial Monte Carlo particle transport study (Geant4-based) shows that in a realistic space mission, the energy range of protons traversing all the three layers (deriving from both the primary space radiation and as products of nuclear fragmentation in the layers) extends until low energies.Low energy protons (up to few hundred keV), which stop inside the cell, have been shown to cause the most severe damage on TJ solar cells [45].For all the reasons listed above, we chose the low energy regime as the focus of this paper.
In this context, the aim of this work is to present a systematic investigation of the ESP of the sub-junction materials Ga 0.5 In 0.5 P, GaAs, and Ge of space solar cells under low-energy (keV range) proton impact in channeling conditions, using RT-TDDFT as implemented in SIESTA [46,33].We consider various channeling trajectories for the projectile (different impact parameters) inside the target.The results are compared with SRIM [5] data, which do not discriminate channelling trajectories from others.
Finally, an analysis is performed of the electronic energy loss at the interface between two upper layers of the TJ structure.To the best of our knowledge, no study has been presented investigating the change of the electronic response to particle radiation across the stack.The results presented here should contribute to understand to what extent channeling effects can be linked to both the electronic density and the atomic number of the element(s) in the compounds constituting the solar cell stack, which should serve as a basis to a further study of non-adiabatic MD simulations of cascades.
The key questions we ask in this work are: (I) how does low-energy channelling electronic stopping changes with respect to the widely used random-trajectory average, as used in SRIM; (II) how large can the differences among channels be; (III) how different are the three materials; (IV) are there significant electronic effects at interfaces; and (V) what is the effect of the strain in the layers due to epitaxial growth.

Method
2.1 Monte Carlo particle transport simulations through the solar cell 1D-Monte Carlo proton transport calculations of the isotropic protons through the multilayer stack of the TJ solar cell have been performed with the MUlti-LAyered Shielding SImulation Software (MULASSIS) [47] (based on Geant4 [8,9,10], available via the European Space Agency's (ESA) online interface SPENVIS [48].Geant4 is a Monte Carlo particle transport simulation software used in application to high energy, nuclear and accelerator physics, as well as studies in medical and space science.In order to simulate different energy ranges, different systems and configurations, a pre-defined set of physics models (the so-called physics lists) is included in the Geant4 toolkit, which allow for the description of hadronic and electromagnetic interactions with different accuracy.The hadron/ion interactions (in general for protons, neutrons, and pions) at energies above 20 GeV, are described by a Quark Gluon String (QGS) model [49].At energies below 10 GeV, the Binary Intra-nuclear cascade (BIC) model for primary proton and neutron interactions with nuclei is used [49].At low and intermediate energies, ionization and elastic scattering are the main processes of the interaction of protons with matter.
The energy loss due to electromagnetic interactions is calculated using the energy loss, range, and inverse range tables pre-computed at initialization of Geant4 for each material [50,51].In all Geant4 electromagnetic physics configurations, two models are used for proton ionisation [52]: the PSTAR/SRIM stopping power [53,5] for proton energies below 2 MeV, and the Bethe-Bloch formula with shell, Barkas and Bloch, and density effect corrections [54] for energies above 2 MeV.The G4EmStandardPhysics − option3 was used for the simulation of electromagnetic interactions.By choosing this physics list, the appropriate descriptions of the electromagnetic processes are automatically applied for different particles in different energy ranges.
A simplified solar cell structure with a typical SiO 2 amorphous 100 µm thick coverglass [55], and with the Ga 0.5 In 0.5 P, GaAs, and Ge layers having 0.8, 2, and 175 µm thickness respectively, is used, similarly to [56].The omnidirectional fluxes of trapped protons, covering the energy range of 0.1 to 400 MeV (Earth's radiation belt, on the basis of data from more than 20 satellites), have been generated via the NASA AP-8 model [57].
The Transport of Ions in Matter (TRIM) part of the SRIM code is used to calculate the proton track trajectories through the solar cell stack.TRIM is a Monte Carlo code that calculates the interactions of energetic ions with amorphous targets using several approximations, such as binary collisions only, an analytical formula for the ion-ion interactions, and the so-called concept of a Free-Flight-Path between the collisions, such that only significant collisions are evaluated [5,6].

Electronic stopping power from RT-TDDFT
The ground state configuration of each system is calculated using the DFT implementation of the SIESTA code (see Appendix A. for details) [46] with the projectile placed at the initial position for each trajectory (specified in section 3.2).
The ground state Kohn-Sham orbitals of the system are obtained by solving the Kohn-Sham equations self-consistently until the total ground state energy is converged.The exchange-correlation functional uses the local-density approximation (LDA) in the Ceperley-Alder form [58]. Norm-conserving Troullier-Martins [59] relativistic pseudopotentials are used to replace the core electrons (see Appendix A. for the parameters used to generate pseudopotentials).It is known that core electrons can have large effects on the ESP (see e.g.Ref. [60]).However, here we only test the effect of core electrons for one proton trajectory in Ge, for all the others we consider the interaction of the projectile with valence electrons, given the low charge of the projectile, and the relatively low velocity range.In such limit the effect of core electrons is negligible as we demonstrate in our work (see also the validation with experiments of similar calculations for protons in Ge in Ref. [33]).The valence electrons are represented by a double-ζ polarized basis set of numerical atomic orbitals defined as specified in the Appendix A.
We use real-time time-dependent DFT (RT-TDDFT) implementation of the SIESTA method [61,33] to evolve the electronic orbitals in time.The time-dependent Kohn-Sham equations are solved by real-time propagation for discretized time, using the Crank-Nicholson scheme [63,61] with a time step dt = 1 attosecond.The effect of the moving basis set is accounted for by a Löwdin transformation, as described in Refs.[26,64], which is known to be adequate for relatively low velocity of the moving basis orbitals, and offers strictly unitary propagation for finite dt.
The focus of this work is on the purely ESP for constant velocity projectiles.The forces on the nuclei of the target atoms and on the projectile itself are therefore disregarded for the time propagation, thereby describing nuclear dynamics with frozen host nuclei and a constant velocity projectile, as done in many similar studies [25,30,33,60].It allows the clean separation of the electronic and nuclear contributions to the total stopping.Fixing the atomic positions has a negligible effect for the simulations performed in this work, given the fact that the projectile traverses the simulation box in about t ≤ 4 fs, and the ionic displacements on that time scale are negligible (head-on collisions are avoided in channelling conditions).
All three target materials have a diamond cubic structure with a similar value of the lattice constant.We used the theoretical lattice constant of 5.59 Å for Ge (for which the experimental value is 5.66 Å) in order to compare our results with Ullah et al. [33].For other layers, we used experimental lattice constants of 5.65 Å f or GaAs, and 5.66 Å for Ga 0.5 In 0.5 P (average between the values for GaP, 5.45 Å, and InP, 5.87 Å).The change in the lattice constant by 1.7% only changes the stopping power by a negligible 0.28% according to our results (see section 3.6).
The targets are modeled by a super-cell of 96 atoms, constructed by multiplying the conventional unit cell (which consists of 8 atoms) by 2 × 2 × 3 in the x, y, and z directions, respectively.A 2 × 2 × 1 Monkhorst-Pack [65] k-point mesh is used, which corresponds to a k-grid inverse cutoff of 8.4 Å.The k-points are displaced to the centers of the grid cells.Two channels in the target super-cell are chosen as trajectories for the proton, along the [001] and [011] crystallographic directions, as shown in figure 2 for the case of Ge.Besides the channeling trajectories, and defining the channel center as the axis furthest away from all atoms, we also choose various off-center parallel channel trajectories with different impact parameters, i.e., at different closest distances of the proton from the host atoms.The same super-cell is used in simulations for [001] and [011] channels, since the length of the trajectory is similar in both directions.
The ESP is calculated as S = dE tot /dz, the derivative of the Kohn-Sham total electronic energy with respect to projectile position z along the constant-velocity path [25].E tot (t) is taken as approximate density-functional for H (t). The average of S is obtained from the slope of a linear fitting to E tot (z) as shown in figure 3. Figure 3 shows the energy for the case of a proton moving in Ge with a velocity of 0.5 a.u.along the center of the [011] channel and for an off-center trajectory along the same direction.The oscillations in the figure reflect the periodicity of the crystal.They are quite prominent for the low-impact factor case, but barely noticeable for the mid-channel trajectory.
All the results in the following are compared with the SRIM [5] data for the ESP.SRIM results are obtained semiempirically by averaging over a number of different incident directions with distinct impact parameters, with no explicit  4 clearly show that the coverglass lowers the fluences of lower energy particles, leaving almost unaffected the high energy portion of the proton spectrum.Those protons are not only the primary impacting protons, but also those generated by the nuclear interactions in the target.Overall, the results show that, for a realistic scenario, once the particles are slowed down by the 100 µm SiO 2 coverglass, the energies of the particles entering all the other layers change to a negligible amount (except for the particles exiting the bottom of the bottom Ge subcell, whose spectrum is considerably "hardened", i.e., energies lower than ∼12 MeV do not exit the whole stack (figure 4).Within the protective amorphous oxide cap layer, ions will be scattered with no channeling, and only a few will find the channels.Contrarily, in the crystalline materials constituting the other layers, channeling conditions may either enhance the depth of penetration or exhibit, for specific impact parameters, blocking trajectories.The fluence always contains particles extending to the low energy regime of one-hundredth to tenths of keV, for essentially all layers in the structure.
Figure 4: Fluence of the omnidirectional primary proton radiation across all layers of the solar cell structure (coverglass included).
In figure 5, we report the results of the SRIM (TRIM) Monte Carlo simulations, which are performed, as it is the case in ground testing studies, assuming unshielded configurations and unidirectional (normal incidence), monoenergetic proton irradiation.TRIM only deals with Coulomb scattering.Thus, the calculated damage is that produced by the projectile itself, the primary knock-on atoms (PKAs), and consequently generated secondary (SKAs), and not by eventual "secondary particles" produced in nuclear reactions.The energies chosen for the impacting protons are representative of the low energy ranges used by recent studies (3 MeV [56], 1 MeV [66, 67, 68], 0.3 MeV [69], and 0.1 MeV [70]).It is important to note that, for such impacting energies, the particles will stop inside the solar cell, in the Ge layer (figure 5a, b), in the GaAs layer (figure 5c), or between the two active layers Ga 0.5 In 0.5 P/GaAs (figure 5d).Thus, the relevant range of impact energies that is necessary to consider in the ESP calculations extends down to 0 keV.
Since the damage produced per unit path length increases as the proton energy decreases and some protons stop within the active region of the cell (figure 5), a non-uniform damage can be induced at such low energies.Thus, the application of the current modeling approach, the so-called DDD (displacement damage dose) model exhibits some limitations for a relevant energy range (1-10 MeV), since it is based on a uniformity of the damage [71] and on an observed linear dependence between the NIEL (i.e., the energy deposited to the atomic displacements) and the damage (i.e., the number of defects produced).Thus the NIEL concept should be improved in order to take into account the details of the structure of the system and the dependence of the NIEL on depth (instead of the NIEL based on the incident energy [70]), which can change in channeling vs non-channeling conditions).If an ion enters a crystalline target along a direction with a low Miller index, the energy-dependent ratio of electronic and nuclear stopping, which is valid for an amorphous material, will no longer be valid.

Electronic stopping power in channeling conditions, comparison to SRIM
The ESP as a function of the proton velocity is shown in figure 6 for three layers of the TJ solar cell and for two channels: [001] and [011].The corresponding incident energies of the proton can be obtained from E = m p v 2 /2, where m p = 1836 a.u. is the proton mass (0.999−20.234 keV for the chosen range of velocities of 0.2−0.9a.u.).The trajectories of the proton in the [001] and [011] channels are schematically depicted in the insets of figure 6. Trajectory 1 corresponds to the center of the channel in all the cases.The impact parameter for the off-center trajectories closest to the rows of atoms (and to the axis between two atoms in the case of the channel [011]) is 0.7 Å (approximately half the distance from the center of the [001] channel to the target atom).In Ge, we also consider additional trajectories with the impact parameter equal to 1.05 Å for trajectory 2 in figure 6a and 1.4 Å and 1.86 Å for trajectories 2 and 4, respectively, in figure 6b.The orange curve in each panel of figure 6 shows the ESP calculated with SRIM [5] using the following values of the density: 5.53 g/cm 3 for Ge, 5.32 g/cm 3 for GaAs, and 5.52 g/cm 3 for Ga 0.5 In 0.5 P, which correspond to the lattice constants used in RT-TDDFT calculations.
In Ge, the values of the ESP are very similar for all three [001] trajectories (figure 6a), since the electronic density distribution does not vary much inside this narrow channel.This can be seen in figure 7 in which we plot the electron density averaged along the [001] channel as a function of position in the perpendicular plane for all three layers.The average is obtained from the density calculated for 160 [x, y] planes along the z-axis with a step of dz = 0.1 Å. Figure 7c shows that the average density is much more homogeneous throughout the unit cell in Ge as compared with the other two layers (figure 7b, c) and that the density varies only slightly along three trajectories.
A larger dispersion of ESP is observed in the [011] channel in Ge (figure 6b), being a consequence of a larger variation of the electron density inside this channel [33].The lowest ESP comes from the mid-channel, as it can be reasonably expected given the lower electron density encountered in such trajectory.For a slow proton moving along the [011] channel, the stopping power is similar for all the trajectories.However, at higher velocities, the stopping power depends on the impact parameter quite significantly.As it has been noted in Ref. [33], this behavior is a consequence of the strong correlation between the stopping power and the average local density within a small radius of the impact parameter, and that this radius is larger for a slower projectile.
On average, our results for the ESP in Ge are slightly lower than the SRIM for the proton moving in both channels, except for the trajectories 3 and 5 in the [011] channel, for which the ESP approaches the SRIM values.This suggests that in the SRIM calculations, the average density is higher than the density on most trajectories in the RT-TDDFT calculations for Ge.Present results for a proton moving on the mid-channel trajectories along the directions [001] and [011] in Ge are in agreement with the results previously obtained by Ullah et al. [33].
In figure 6a, we show as well the ESP for an off-channel proton trajectory, which makes 14 • with the direction [001] (z−axis).The initial projectile position is at the center of the [001] channel at a distance 0 Å from the first atomic plane perpendicular to this direction.The velocity of the proton is 0.62 a.u.(0.6 a.u.along the z−axis and 0.15 a.u.along the x−axis).The calculation was performed during 6 fs, four times the time needed to cross the simulation box at this The colour scale is the same for all three panels.
velocity.Thus it is equivalent to four different random trajectories [32].The resulting ESP for the off-channel trajectory is in between the values for two off-center channeling trajectories 2 and 3, and thus does not bring the channeling results closer to SRIM.Possible reasons for that is the absence of core electrons, as discussed in the next subsection, and insufficient trajectory sampling when calculating the ESP for the off-channelling case (i.e., a longer propagation is probably needed).Other reasons may come from the inaccuracy of SRIM.
In GaAs, similarly to Ge, the ESP for the proton in the [001] channel is almost independent on the impact parameter and agrees perfectly with the SRIM data (figure 6c).Only the ESP for the proton on trajectory 3 is slightly higher then on 1 and 2, given the higher average electron density on this trajectory (figure 7b).In the [011] channel, the ESP is below SRIM for the proton moving on the mid-channel trajectory, while it is above SRIM for all the off-center trajectories with small impact parameter (i.e., higher electron density) (figure 6d).At proton velocities below 0.5 a.u., the ESP in the [011] channel is similar for all the off-center trajectories.At higher velocities, however, the ESP corresponding to the trajectory 3, closest to the As atoms, is slightly higher.This could indicate that the maximum of the ESP is located at different velocities for different proton trajectories.
In the case of the Ga 0.5 In 0.5 P layer, again, on average, our results are in good agreement with SRIM (figure 6e, f ).However, the dependence of the ESP on the trajectory is observed in both channels.For the proton moving on trajectories close to Ga and In atoms in both channels, the ESP values are very close, owing to the fact that both elements have three electrons in the outer shell.A much higher ESP is obtained for the trajectories close to P, which have five valence electrons, and thus, the electron density is higher along these trajectories (trajectory 3 in figure 7a).Analysis of the average electron density for different proton trajectories in the direction [001] in Ga 0.5 In 0.5 P, indicated with open circles in figure 7a, shows that the density is similar for trajectories 1, 2, and 4, but is higher by a factor of 2 along the trajectory 3.However, the stopping power for a proton moving on the trajectory 4 is only 1.4 (maximum, at v = 0.7 a.u.) times larger than the one on the trajectory 1 (figure 6e).Thus, the stopping power is not linear with the average density along the proton path.
In figure 6e, f we also compare our results to the calculations by Lee et al. [32].Similarly to our work, the ESP is obtained with RT-TDDFT (LDA), but using a plane-wave basis set.They are represented as empty squares for the stopping along the mid-channel trajectory, and as empty diamonds for the off-channel trajectory (i.e., random trajectory through the lattice, not parallel to any of the main crystallographic directions).There are small deviations between the results of both approaches, probably due to basis sets and pseudopotentials, but the main results, trends and conclusions are the same.Those deviations are smaller than the spread for different trajectories, and the differences with respect to SRIM.
Small deviation of the ESP from the SRIM behavior observed at high velocities (especially in Ge) is due to two reasons.First, closer to the maximum of the stopping power, the effect of the inner electrons starts to be important, which is not taken into account in these calculations.We will discuss this effect in the next subsection.Another reason is that the Sankey integrator used in SIESTA to propagate the KS states in time [64] is known to be problematic at higher energies.

Comparative analysis of the ESP in different layers
We further analyze our results by comparing the values of the ESP for equivalent proton trajectories inside different layers of the TJ solar cell as shown in figure 8. Overall, the ESP do not vary much from material to material for the proton moving on the mid-center trajectories (figure 8a, b).For the smallest impact parameter of 0.7 Å in the [001] channel (figure 8c), we show the ESP for a proton trajectories close to different elements in each material.The highest ESP corresponds to the trajectory close to P atoms in Ga 0.5 In 0.5 P. Figure 8d shows the ESP for the proton trajectories 5 (Ge) and 4 (GaAs and Ga 0.5 In 0.5 P) in the [011] channel (see trajectories in figure 6).Here as well, the largest ESP is obtained for the proton on the trajectory close to P atoms.
The importance of core electrons was checked for the case of Ge (figure 8a) by including 3d electrons into the valence shell (see Appendix for details of the calculations).At velocities below 0.6 a.u., the addition of core electrons has no effect on the ESP.However, at velocity 0.8 a.u, the ESP is notably higher when 3d electrons are included in the valence, for this particular trajectory.
Since Ga 0.5 In 0.5 P is a compound with the largest number of constituent species in TJ solar cells, it is interesting to look at the electron density distribution inside it and try to correlate it with the observed ESP. Figure 9 shows the electron density for three different planes perpendicular to the [001] direction.Each plane crosses the positions of atoms of only one type (Ga, In, or P).From figure 9, it is obvious that the electron density in the vicinity of the P atoms is much higher, which explains the largest ESP for protons on the trajectories 3 and 4 (figure 6).
Not only the value, but the spatial distribution of the electron density varies in different planes as well.Figure 10 shows the density as a function of y coordinate for all three panels of figure 9 corresponding to x = 0. Thus three curves in figure 10 show the density distribution around one single atom of Ga, In, and P, respectively.The system of coordinates is chosen so that all the atoms are centered at the origin (y = 0).The density, going down towards the nuclear position, reflects the absence of core electrons in the calculations.The maximum density around P atom is higher and is closer to the origin.Thus the projectile moving at 0.7 Å from the atomic position (indicated by the vertical dashed line in figure 10), in the case of P, is moving in the area of much higher electron density than in the case of Ga and In.

Comparison with the homogeneous electron gas (HEG) model
In order to show the non-linearity of the ESP obtained with RT-TDDFT, we compare our results with the HEG approximation.In HEG, the electron density is defined through the so-called Wigner-Seitz radius r s , a one-electron radius, as n = 3/(4πr 3 s ).The stopping power in HEG is calculated as a product of the friction coefficient and proton velocity.The friction coefficients corresponding to different r s are taken from [72].In figure 11 we show the stopping power as a function of the average electron density in Ga 0.5 In 0.5 P channel [001].The lower density points (three points close to each other) correspond to the trajectories 1, 2, and 3, and the higher density corresponds to the trajectory 4.
The crosses show the results of the HEG approximation, while the squares show the RT-TDDFT stopping power.The comparison is shown for three different velocities of the proton.At lower velocities, the agreement is very good.At higher velocity, however, the importance of the non-linear effects is clearly observed as the HEG results overestimate our results.

Interface effects on the proton energy loss
In order to study the effect of an interface on the energy loss of the proton, we constructed a super-cell of the interface between the two upper layers of the TJ solar cell, i.e., Ga 0.5 In 0.5 P/GaAs along the crystallographic direction [001] of the epitaxial growth of the lattice-matched solar cells (figure 12a).The super-cell consists of 192 atoms and has the lattice parameter average between the pure bulk Ga 0.5 In 0.5 P and GaAs.We have chosen two trajectories of the proton through the interface, the channeling one (figure 12b) and the off-center channeling trajectory (figure 12c), in which the proton first passes close to P and then to As atomic rows in the direction [001] with an impact parameter of 0.7 Å.
The difference in the energy loss of the proton moving through the interface and through the pure Ga 0.5 In 0.5 P and GaAs super-cells of the same size is shown in figure 13 for velocities 0.3 a.u. and 0.6 a.u.for both trajectories.The interface between two materials is located at z int = 16.26Å.The proton first moves though the Ga 0.5 In 0.5 P, hence the difference in the proton energy loss ∆(dE tot (z)) = dE int (z) − dE GaInP2 (z) is equal to zero until z z int in all the cases.For a similar reason, starting from z z int , the value of Interestingly, the slope of ∆(dE tot (z)) has opposite trends for the two trajectories.In the case of the mid-channel trajectory (figure 13a, c), the proton energy loss in GaAs is larger than in Ga 0.5 In 0.5 P. On the off-center trajectory (figure 13b, d), however, the behavior is reversed.This change of slope occurs due to the fact that the electron density around P atoms is more localized (see figure 9), and is slightly lower than around Ga in the center of the [001] channel (y = 2.1 Å, figure 10).Thus, the electron density in the center of the channel in Ga 0.5 In 0.5 P is expected to be lower than in GaAs, giving rise to lower energy loss of the proton moving on this trajectory in Ga 0.5 In 0.5 P. On the off-center trajectory, on the contrary, the density around P is much higher.This behaviour is, however, a bulk effect not related to the interface itself.
To quantify the interface effect and separate it from the bulk effects, we have calculated the proton energy loss on a given part of the trajectory equal to the length of two unit cells on each side of the interface structure (dE int(Ga0.5In0.5P)and dE int(GaAs) ) and for the equivalent part of the trajectory of pure structures (dE Ga0.5In0.5Pand dE GaAs ).The difference between the energy loss of a proton moving through the interface structure and the sum of the energy losses of a proton moving through the corresponding intervals of pure structures, gives us the difference due to the interface: ∆E int = dE int(Ga0.5In0.5P)+ dE int(GaAs) − dE Ga0.5In0.5P− dE GaAs .The values we obtained for ∆E int for center and off-center channeling trajectories are 0.39 eV and 0.60 eV for velocity 0.3 a.u., and −0.075 eV and 0.072 eV for velocity 0.6 a.u., respectively.Thus, the interface effect is higher at lower velocities.The energy loss due to the interface even becomes negative at high velocity in the case of center-channel trajectory, meaning effectively an energy gain.These differences, however, are 1 − 2 orders of magnitude smaller than the difference between the materials, i.e., the bulk effect.

Effect of strain in lattice-matched solar cells on the ESP
Most of the TJ solar cells currently used for space applications are lattice-matched.Ge is the thickest layer used as a substrate.The GaAs layer is grown on top of Ge and the Ga 0.5 In 0.5 P is grown on top of the GaAs layer, with the lattices of different layers matching at the interface.The experimental lattice constants of these materials are very similar, 5.658 Å for Ge, 5.653 Å for GaAs (only 0.08% different from Ge), and 5.6596 Å for Ga 0.5 In 0.5 P.This small difference in the lattice constant leads to a strain in the GaAs and Ga 0.5 In 0.5 P layers.We have analysed the effect of strain in two cases.In the first case, we have calculated the stopping power for a proton moving with v = 0.8 a.u in GaAs with lattice constants 5.65 and 5.75 Å (the difference of 1.74%).The stopping power only changed by 0.28%.In another case, we compared the stopping power for a proton moving with v = 0.5 a.u in GaAs with lattice constants 5.75 and 5.57 Å (3.13% different).The difference in the stopping power was only 1.13%.The actual differences in the lattice constants of the three layers of TJ solar cell are much smaller, and the changes in the stopping power can be considered negligible.

Conclusions
In this work, using RT-TDDFT we have calculated the electronic stopping power of the three sub-junctions of the TJ solar cells for the impacting protons.Kinetic energies of the proton are considered in the keV range, which we have shown to be relevant in the transport of protons through the full stack in a realistic scenario.We compared our results for channeling conditions to the SRIM semi-empirical stopping power.Our results have shown that in Ge, the stopping power in the channel [001] is almost independent on the trajectory, while in both compounds (GaAs and Ga 0.5 In 0.5 P) the electronic stopping varies more along different trajectories, affected by varying electron density.In the case of the channel [011], strong dependence of the stopping power on the trajectory is observed in all three materials.The effect is more significant at higher proton velocities.
The change of sign for the slope of the energy loss difference between GaAs and Ga 0.5 In 0.5 P in the interface structure is a bulk effect, which is explained in terms of electron density.Namely, the average density along the center of the [001] channel in GaAs is higher than in Ga 0.5 In 0.5 P, while it is the other way around along the off-center-channel trajectory.
The effect of the interface between the layers of the lattice-matched multilayer solar cell, as well as the effect of strain, have been found to be negligible, which is understood, given the very similar chemistry and lattice constants of the three materials.The interface effect is found to be higher for lower velocity.At higher velocity, the energy loss can become negative, which effectively means that the proton loses less energy due to the interface.
Additional studies can help to understand the importance of the channeling effect in reducing the damage caused by radiation.The coupling between nuclear and electronic stopping power will be considered as a further step in this study.

Figure 2 :Figure 3 :
Figure 2: Ge super-cell (96 atoms).A proton position is shown in the center of (a) [001] channel and (b) [011] channel.GaAs and Ga 0.5 In 0.5 P have a similar structure.

Figure 6 :
Figure 6: Electronic stopping power for a proton in Ga 0.5 In 0.5 P/GaAs/Ge compared to SRIM.The stopping power is shown as a function of velocity for the channel [001] in (a) Ge, (c) GaAs, and (e) Ga 0.5 In 0.5 P; and for the channel [011] in (b) Ge, (d) GaAs, and (f) Ga 0.5 In 0.5 P. SRIM results are shown as the solid curve with crosses for each layer.A schematic representation of the trajectories of the proton in each channel are shown in the insets.Green asterisk in (a) shows our result for an off-channel trajectory (see text for details).Empty symbols in (e) and (f) are the results from Lee et al. [32], the squares and diamonds correspond to the channeling and off-channeling results, respectively (see text for details).

Figure 7 :
Figure 7: Electron density averaged along the [001] channel as a function of position in the perpendicular plane for Ga 0.5 In 0.5 P, GaAs, and Ge.The trajectories of the proton in all systems are indicated with open circles and are the same as in the insets of figure 6.(a) average density in Ga 0.5 In 0.5 P, black dots indicate positions of the alternating Ga/In atoms, while black stars indicate positions of the P atoms; (b) average density in GaAs, black dots and squares show positions of the Ga and As atoms, respectively; (c) average density in Ge, black dots show positions of the Ge atoms.The colour scale is the same for all three panels.

Figure 8 :
Figure 8: Electronic stopping power for a proton in Ga 0.5 In 0.5 P/GaAs/Ge, a comparison of different materials.The stopping power is shown as a function of velocity for (a) direction [001], center of the channel.Red diamonds show the results for Ge including semicore electrons; (b) direction [011], center of the channel; (c) direction [001], impact parameter of 0.7 Å; (d) direction [011], impact parameter of 0.7 Å from the edge of the hexagon (see insets of Fig.6 (b,d,f)).

Figure 9 :
Figure 9: Electron density for three planes in Ga 0.5 In 0.5 P perpendicular to the direction [001].The planes are crossing through (a) Ga; (b) In; (c) P atoms.

Figure 10 :
Figure 10: Electron density along the y-axis, corresponding to the cut through x = 0 in three panels of figure 9.The vertical dashed line shows the position of the proton in the off-center channeling trajectory at 0.7 Å from the atomic position, corresponding to trajectories 2, 3, and 4 in the inset of figure 6e.

Figure 12 :
Figure 12: Atomic structure of the interface.(a) Structure of the interface between Ga 0.5 In 0.5 P and GaAs, (b) proton in the center of the channel [001], (c) proton in the off-center channel at 0.7 Å from P and As atoms of the interface structure.

Figure 13 :
Figure 13: Difference in the proton energy loss between pure GaAs and Ga 0.5 In 0.5 P and the interface between the two materials.The quantity plotted is ∆(dE tot (z)) = dE int (z) − dE pure (z), where dE(z) = E(z) − E(z = 0) for (a) channeling trajectory, v = 0.3 a.u.; (b) impact parameter 0.7 Å, v = 0.3 a.u.; (c) channeling trajectory, v = 0.6 a.u.;(d) impact parameter 0.7 Å, v = 0.6 a.u., as indicated in each panel.The interface between GaAs and Ga 0.5 In 0.5 P is located at 16.26 Å and is showed by the vertical dashed line in each panel.