Dimensionality and Anisotropicity Dependence of Radiative Recombination in Nanostructured Phosphorene

TL;DR

This paper develops a first-principles framework to study how dimensionality and anisotropicity influence radiative recombination in nanostructured phosphorene, revealing significant effects on exciton lifetime and polarization of emitted light.

Contribution

It introduces a general first-principles approach to analyze exciton dynamics across different dimensions in low-dimensional semiconductors, highlighting the impact of anisotropicity.

Findings

Radiative lifetime increases by about 1000 times when reducing dimensionality from 2D to 0D.

Finite temperature radiative lifetime calculations require accurate exciton dispersion models.

Monolayer phosphorene emits linearly polarized light, unlike isotropic 2D materials.

Abstract

The interplay between dimensionality and anisotropicity leads to intriguing optoelectronic properties and exciton dynamics in low dimensional semiconductors. In this study we use nanostructured phosphorene as a prototypical example to unfold such complex physics and develop a general first-principles framework to study exciton dynamics in low dimensional systems. Specifically we derived the radiative lifetime and light emission intensity from 2D to 0D systems based on many-body perturbation theory, and investigated the dimensionality and anisotropicity effects on radiative recombination lifetime both at 0 K and finite temperature, as well as polarization and angle dependence of emitted light. We show that the radiative lifetime at 0 K increases by an order of $10^3$ with the lowering of one dimension (i.e. from 2D to 1D nanoribbons or from 1D to 0D quantum dots). We also show that obtaining the radiative lifetime at finite temperature requires accurate exciton dispersion beyond the effective mass approximation. Finally, we demonstrate that monolayer phosphorene and its nanostructures always emit linearly polarized light consistent with experimental observations, different from in-plane isotropic 2D materials like MoS2 and h-BN that can emit light with arbitrary polarization, which may have important implications for quantum information applications.

Tables

Table 1. Table 1: Radiative lifetime of the first exciton of MBP with 2D, 1D and 0D dimensions at 0 K and 263 K, with light polarized along the armchair direction.

System	Width Armchair (Å)	Width Zigzag (Å)	$\hslash {\Omega }_0$ (eV)	${\mu }_A^2$ (a.u.)	Lifetime at 0 K
0D-2a2z	9	9	2.80	0.75	57 ns
0D-3a2z	13	9	2.39	3.49	20 ns
0D-4a2z	18	9	2.16	6.43	14 ns
1D-a-2z	$\infty$	9	1.73	7.42	21 ps
1D-a-3z	$\infty$	12	1.52	9.78	21 ps
1D-a-4z	$\infty$	15	1.51	11.89	17 ps
1D-z-2a	9	$\infty$	2.72	0.04	1.0 ns
1D-z-3a	13	$\infty$	2.40	0.45	130 ps
1D-z-4a	18	$\infty$	2.17	0.64	110 ps
2D*	$\infty$	$\infty$	1.54	2.56	99 fs
2D(Exp)*					211 psYang et al. (2015)
2D(noEM)*					101 ps
2D(EM)*					1.8 ns

Equations

$$\gamma (Q_{\text{ex}},q_{\text{L}},\lambda )$$

$$\times \delta (E(Q_{\text{ex}})-\hslash c{q}_L),$$

$$\gamma (Q_{\text{ex}})$$

$$\text{matrixelement}G,1_{q_{\text{L}}\lambda }{H}^{\text{int}}S(Q_{\text{ex}}),0$$

$$=$$

$$\gamma (Q_{\text{ex}})$$

$${\Omega }_{0}l/c\approx 10^{-2}-10^{-3}.$$

$$\gamma (T)$$

$$Z$$

$$E_i(Q_{\text{ex},i})=$$

$$+\frac{4\pi }{d}{\mu }_i^2\left(1-\frac{1-e^{-Q_{\text{ex},i}d}}{Q_{\text{ex},i}d}\right),$$

$$E(Q_{\text{ex}})=$$

$$+(E_y(Q_{\text{ex},y})-E_0).$$

$$I(q_{\text{L}},\lambda )$$

$$I({\theta }_{\text{L}},{\phi }_{\text{L}},\lambda )$$

$$I({\theta }_{\text{L}},{\phi }_{\text{L}},\text{IP})$$

$$I({\theta }_{\text{L}},{\phi }_{\text{L}},\text{OOP})$$

$$-M_y\cos {\theta }_{\text{L}}\sin {\phi }_{\text{L}}+M_z\sin {\theta }_{\text{L}}{)}^2,$$

$$S_0$$

$$S_1$$

$$S_2$$

$$S_3$$

$$S_3=n_0{\Gamma }_{0}i\cos {\theta }_{\text{L}}(M_x^{\ast }{M}_y-M_x{M}_y^{\ast }).$$

Full text

Two-dimensional (2D) materials are known for their unprecedented potential in ultrathin electronics, photonics, spintronics and valleytronics applicationsLiu et al. (2017); Youngblood et al. (2015); Shiue et al. (2017); Guo et al. (2016); Bernardi et al. (2017); Eda and Maier (2013); Britnell et al. (2013); Wong et al. (2017); Han (2016); Feng et al. (2017); Avsar et al. (2017). Strong light-matter interactions and atomically-thin thickness lead to exotic physical properties which do not exist in traditional 3D semiconductors. In particular, the monolayer black phosphorus (MBP), commonly known as phosphorene, has recently attracted significant interest due to its emerging optoelectronic applications and the development of large scale fabrication methodsCupo and Meunier (2017); Akhtar et al. (2017); Carvalho et al. (2016).

Unlike other common 2D materials such as graphene, $h$-BN, and transition metal dichalcogenides, MBP has a strong in-plane anisotropic behavior along two distinctive directions, denoted as “armchair”and “zigzag” Cupo and Meunier (2017); Pang et al. (2017); Li et al. (2017); Tran et al. (2014); Qiu et al. (2017). For example, the lowest exciton transition is only bright when light is polarized along the armchair direction and the corresponding excitonic wavefunction is much more extended along this direction, leading to so-called “quasi-1D” excitonsTran et al. (2014). The complex nature of excited states and the anisotropic behavior of excitonic wavefunctions may lead to intriguing effects on the exciton dynamics in MBP, which are distinctive from the in-plane isotropic 2D systems as we will discuss later. Further lowering dimensionality is expected to show an interplay between quantum confinement and quasi-1D excitonic nature which affects its optoelectronic properties in a complex manner. However, because lower dimensional MBP nanostructures (such as nanoribbons and quantum dots) are more difficult to stabilize than two-dimensional MBPCastellanos-Gomez et al. (2014); Yang et al. (2015), the optical measurements and determination of excited state lifetime of these nanostructures are more challenging. Only the exciton recombination lifetime of 2D MBP has been determined to be 211 ps through time-resolved photoluminescence measurements with light polarized along the armchair direction at -$10^{\circ}$CYang et al. (2015). This calls for theoretical studies to predict exciton dynamics of low dimensional MBP nanostructures.

Previously, the radiative lifetime of excitons of typical in-plane isotropic 2D systems have been computed by coupling Fermi’s golden rule with model Hamiltonians or the Bethe-Salpeter equation Palummo et al. (2015); Chen et al. (2018); San-Jose et al. (2016); Ayari et al. (2018). However, the radiative lifetime and light emission intensity of in-plane anisotropic systems (e.g. phosphorene), which have unique dependence on the polarization direction of light, have not been investigated in-depth. Additional outstanding questions involve the dimensionality (i.e. going from the 2D MBP to 1D and 0D nanostructures) and temperature dependence of excited state lifetime. In particular, temperature effects are determined by the exciton dispersion in momentum space, where the typically used effective mass approximation may break down for low dimensional systems Cudazzo et al. (2016). By answering the above questions, we will propose general principles and pathways of engineering radiative lifetime in anisotropic low-dimensional systems.

The rate of emission of a photon with specific wave-vector ${\bf\it q}_{\textrm{L}}$ and polarization direction $\lambda$ from an exciton with wave-vector ${\bf\it Q}_{\textrm{ex}}$ is

[TABLE]

where [math] and $1_{\lambda{\bf\it q}_{\textrm{L}}}$ denote absence and presence of a photon, respectively; $G$ is the ground-state; $S({\bf\it Q}_{\textrm{ex}})$ is the exciton state; $E({\bf\it Q}_{\textrm{ex}})$ is the exciton energy and $V$ is the system volume.

Two important quantities can be derived from $\gamma({\bf\it Q}_{\textrm{ex}},{\bf\it q}_{\textrm{L}},\lambda)$. One is the radiative decay rate $\gamma({\bf\it Q}_{\textrm{ex}})$ of the exciton that only depends on exciton wavevector ${\bf\it Q}_{\textrm{ex}}$; the other is the light emission intensity $I({\bf\it q}_{\textrm{L}},\lambda)$, which only depends on photon wavevector ${\bf\it q}_{\textrm{L}}$ and polarization $\lambda$ for a given system. Both quantities can be directly measured from experiments.

The radiative decay rate of a specific exciton with wave-vector ${\bf\it Q}_{\textrm{ex}}$ can be obtained from:

[TABLE]

the corresponding lifetime can be simply computed as the inverse of $\gamma({\bf\it Q}_{\textrm{ex}})$.

By using the dipole approximation and the relation ${\bf\it p}=-m\frac{i}{\hbar}[{\bf\it r},H]$, the exciton transition matrix element in Eq. 1 can be written explicitly as

[TABLE]

The radiative decay rate $\gamma({\bf\it Q}_{\textrm{ex}})$ can be obtained by combining Eqs. 1-3. The exciton energy $E({\bf\it Q}_{\textrm{ex}})$ and exciton dipole moments $\matrixelement{G}{{\bf\it r}}{S}$ necessary as inputs are computed by solving the Bethe-Salpeter equationPing et al. (2013), which accurately takes into account the electron-hole interaction (this is mandatory for 2D semiconductors where large exciton binding energy $>$ 0.5 eV has been observedHüser et al. (2013); Qiu et al. (2016)).

In order to study the dimensionality dependence of radiative lifetime in MBP nanostructures, we derived the general radiative decay rate expressions for anisotropic materials from 2D to 0D; detailed derivations can be found in Supporting Information (SI). In general, the radiative decay rate can be written as

[TABLE]

where $\gamma_{0}$, which does not depend on the directions of photon or exciton wavevectors, is the exciton decay rate at ${\bf\it Q}_{\textrm{ex}}=0$; the dependence on the wavevector direction is instead contained in $Y({\bf\it Q}_{\textrm{ex}})$, which satisfies $Y({\bf\it Q}_{\textrm{ex}}=0)=1$.

To understand the effects of anisotropicity and dimensionality of nanostructured MBP, we computed electronic structure, absorption spectra and radiative lifetimes of various MBP nanostructures from 2D to 0D as shown in Fig. 1 and Table 1. We used open source plane-wave code Quantum-EspressoGiannozzi et al. (2009, 2017) with Perdue-Burke-Ernzehorf (PBE) exchange-correlation functional Perdew et al. (1996), ONCV norm-conserving pseudopotentials Hamann (2013); Schlipf and Gygi (2015). The band structure is computed from GW approximation with the WEST-code  Govoni and Galli (2015); Pham et al. (2013). The absorption spectra and exciton properties are computed by solving the Bethe-Salpeter equation (BSE) in the Yambo-codeMarini et al. (2009) for MBP nanoribbons and monolayer systems. For MBP quantum dots, a BSE implementation without explicit empty states and inversion of dielectric matrix is applied insteadRocca et al. (2010); Ping et al. (2012); Rocca et al. (2012) to speed up the convergence. More computational details can be found in SI.

Optical absorption spectra of the nanostructures considered in this work are provided in SI; here we will summarize the main features that are relevant to the discussion of lifetimes. The first exciton in phosphorene is bright with light polarized along the armchair direction, but dark along the zigzag direction. This behavior is inherited by phosphorene nanostructures (nanoribbons and quantum dots). However, for these systems the light absorption and dipole matrix elements along a non-periodic direction are significantly weakened by a “depolarization effect”, which comes from the microscopic electric fields induced by polarization charge in an external field. This effect is included in our absorption spectra calculations by taking into account the local field effects in the Bethe-Salpeter EquationMarinopoulos et al. (2003); Bruno et al. (2005).

The radiative lifetimes for various MBP nanostructures at ${\bf\it Q}_{\textrm{ex}}=0$ (obtained by $1/\gamma_{0}$) are summarized in Table 1. The lifetime of the first exciton at 0 K changes dramatically from 0D ($\approx$ 10-100 ns) to 1D ($\approx$ 10-100 ps) then to 2D ($\approx$ 100 fs), and generally decreases with increasing system sizes. Comparing the radiative decay rate expressions in SI we can see that with every decrease of dimension (e.g. from 2D to 1D), an additional $\Omega_{0}l_{i}/c$ multiplicative factor appears in the decay rate expression, where $l_{i}$ is the unit cell size and $\Omega_{0}$ is the exciton energy at ${\bf\it Q}_{\textrm{ex}}=0$. This is due to the momentum conservation between photon $q_{i}$ and exciton $Q_{i}$ along the periodic direction $i$ and energy conservation ($\hbar cq_{L}=\hbar\Omega_{0}$).

Considering typical values in a semiconductor solid such as $\hbar\Omega_{0}\approx 0.5-5$ eV and $l=5-10$ Bohr unit cell lattice parameter along one dimension, we can estimate that

[TABLE]

Accordingly, each reduction of dimension will give approximately $10^{3}$ times longer lifetime at 0 K (this qualitative estimation does not take into account the change of dipole matrix elements and exciton energy). The above discussion explains the increase of the lifetime by several orders of magnitude when decreasing the dimensionality.

For systems of the same dimensionality, the lifetime also varies with the size and periodic direction of nanostructures. Specifically, the radiative lifetime is inversely proportional to the product of squared dipole matrix element ($\mu_{A}^{2}=|\matrixelement{G}{{\bf\it r}}{S}|^{2}$) and the exciton energy ($\Omega_{0}$), i.e. $1/\gamma_{0}\propto 1/(\Omega_{0}^{n}\mu_{A}^{2})$, where $n$ depends on the dimensionality. We found that in most cases, by changing the system sizes, the change of $\mu_{A}^{2}$ is much larger than that of $\Omega_{0}$ and dominates the lifetime, as shown in Table 1. For example, the fact that the radiative lifetime of nanoribbons periodic along zigzag (“1D-z”) is $10^{1}-10^{2}$ times longer than their counterparts periodic along armchair (“1D-a”) is directly determined by $\mu_{A}^{2}$. Because of the anisotropicity of MBP structure, the only non-zero component of $\matrixelement{G}{{\bf\it r}}{S}$ of the first exciton is along the armchair direction. In 1D-a nanoribbons the armchair direction is periodic, so the quasi-1D exciton $\mu_{A}^{2}$ is similar to 2D. On the other hand, in 1D-z nanoribbons the armchair direction is not periodic and $\mu_{A}^{2}$ is strongly reduced due to the depolarization effect; however, the dipole moments with light polarized along the armchair direction is still the main contribution to $\mu_{A}^{2}$. This is different from previous work on 1D nanostructures such as carbon nanotubesSpataru et al. (2005) where only the component $\mu_{z}$ in the periodic direction dominates $\mu_{A}^{2}$. Furthermore, the quasi-1D nature of the first exciton (extended along the armchair direction and confined along the zigzag in Fig. 1) introduces a much larger change in lifetimes as a function of the width for 1D-z nanoribbons compared to 1D-a nanoribbons.

The radiative lifetimes of 0D systems are significantly longer with respect to the other nanostructures, ranging from $2.0\times 10^{4}$ ps to $1.2\times 10^{5}$ ps due to the dimensionality dependent factor for decay rates in Eq. 5 (see Table 1). Similarly to the 1D-z nanoribbons, the 0D-$N$a2z ($N$=2,3,4) series shows a clear trend of lifetime decreasing by increasing the quantum dot size along the armchair direction.

Next we will discuss the exciton recombination lifetime at finite temperature for 2D phosphorene in order to compare with the experimental results. The radiative lifetime at finite temperature $T$ is computed by assuming that the recombination process is slow enough to allow excitons at different ${\bf\it Q}_{\textrm{ex}}$ to reach thermal equilibrium Spataru et al. (2005). The temperature-dependent radiative lifetime is then written as a thermal statistical average: Spataru et al. (2005); Palummo et al. (2015)

[TABLE]

where $E_{0}=E({\bf\it Q}_{\textrm{ex}}=0)$, and $Z$ is the partition function and the condition $Q_{\textrm{ex}}<E_{0}/\hbar c$ is imposed by energy conservation. The evaluation of Eqs. 6-7 requires the knowledge of the exciton energy $E$ as a function of ${\bf\it Q}_{\textrm{ex}}$.

In previous work Spataru et al. (2005); Palummo et al. (2015); Chen et al. (2018) the effective mass approximation for exciton dispersion $E({\bf\it Q}_{\textrm{ex}})=E({\bf\it Q}_{\textrm{ex}}=0)+\hbar^{2}Q_{\textrm{ex}}^{2}/2m$ has been used. However, the effective mass approximation works best for Wannier-Mott excitons in 3D bulk material. In low dimensional materials such as phosphorene, the weak dielectric screening and delocalized wavefunctions in plane introduce strong long-range electron-hole exchange interactions. Thus the exciton dispersion violates the parabolic relation for small values of $Q_{\textrm{ex}}$ Cudazzo et al. (2015, 2016); Koskelo et al. (2017). This calls for a more sophisticated approach to deal with the thermal distribution of excitons beyond the effective mass approximation.

Following Ref. 25, we fit $E({\bf\it Q}_{\textrm{ex}})$ along the armchair and zigzag directions of phosphorene independently using an equation based on the combination of localized and delocalized exciton models:

[TABLE]

where $i$ denotes the direction, $m$ is the effective mass, $\mu_{i}$ is the component of dipole moment matrix along the specific direction $i$, and $d$ is the thickness of the assumed dielectric medium. The partition function is then computed numerically from the fitted $E({\bf\it Q}_{\textrm{ex}})$ assuming that the dispersion along two directions is not coupled:

[TABLE]

The temperature-dependent radiative lifetimes of phosphorene (specifically at 263 K and 0 K) are listed in Table 1. If the traditional approach based on the effective mass (“EM” value in Table 1) approximation is used, a lifetime of 1.8 ns is obtained, which overestimates the experimental result by an order of magnitude. If the more complex model in Eq. 8 is used (“noEM” value in Table  1), the computed lifetime of phosphorene at 263 K is 123 ps, which is in good agreement with the experimental result 211 psYang et al. (2015). The remaining discrepancy (less than a factor of 2) could be due to the use of a Si/\ceSiO2 substrate in experimental measurements, which introduces an additional dielectric screening in the material compared to our free-standing system. Our findings show the importance of the accurate exciton dispersion $E({\bf\it Q}_{\textrm{ex}})$ beyond the effective mass approximation.

After the discussions of radiative decay rates of excitons that provide time-resolved information, we are going to investigate the polarization-resolved and angle-resolved information of the light emitted from the exciton recombinations. Angle-resolved and polarization-resolved light emission measurements are performed with a detector that collects photon with momentum ${\bf\it q}_{\textrm{L}}$ and polarization $\lambda$ within a small area around a spherical angle. The intensity $I({\bf\it q}_{\textrm{L}},\lambda)$ emitted from all possible excitons is

[TABLE]

where $n({\bf\it Q}_{\textrm{ex}})$ is the occupation number of the exciton state at momentum ${\bf\it Q}_{\textrm{ex}}$.

By integrating over the norm of ${\bf\it q}_{\textrm{L}}$ and keeping only the non-zero $\gamma({\bf\it Q}_{\textrm{ex}},{\bf\it q}_{\textrm{L}},\lambda$) in the summation of ${\bf\it Q}_{\textrm{ex}}$ in Eq. 10 (based on the momentum conservation between ${\bf\it Q}_{\textrm{ex}}$ and ${\bf\it q}_{\textrm{L}}$ along the systems’ periodic directions), the light emission intensity along a specific direction of photon wavevector ($\theta_{\textrm{L}},\varphi_{\textrm{L}}$; see Fig.1 in SI) and polarization ($\lambda$) can be expressed as follows:

[TABLE]

where ${\bf\it M}$ is the normalized exciton dipole moment.

We consider intensity along two polarization directions:

[TABLE]

where IP (in-plane) and OOP (out-of-plane) denote two polarization directions that are both perpendicular to the photon wavevector ${\bf\it q}_{\textrm{L}}$ (see Fig.1 in SI), and $M_{i}=\mu_{i}/\mu_{A}$ representing the components of the normalized dipole moment along the directions $i=x,y,z$. The formulation in Eqs. 11-13 is general and only the angle independent term $\Gamma_{0}$ describes the difference among 2D, 1D and 0D nanostructures (detailed derivations can be found in SI).

The emitted photon polarization can be described by Stokes parametersHecht (2002); McMaster (1954, 1961); Peach et al. (2009) $(S_{0},S_{1},S_{2},S_{3})$, which can be written explicitly as the following:

[TABLE]

where \lambda=$$$ bisects the angle between IP and OOP, and \lambda=L$denotes the left circularly polarized light. When$|S_{3}/S_{0}|=0(1), the light is fully linearly polarized (circularly polarized). By evaluating Eq. [14](#S0.E14) with inputs from Eqs. [11](#S0.E11)-[13](#S0.E13), we can get the Stokes parameters as a function of normalized exciton dipole moments{\bf\it M}$.

We found that because of the in-plane anisotropicity in MBP and its nanostructures, the emitted light is completely linearly polarized. As we can consider $M_{z}\approx 0$ (where z is along the perpendicular direction to the material plane), the circular polarization component $S_{3}$ in Eq. 14 reduces to the following form (detailed derivations can be found in SI):

[TABLE]

When one of the following conditions is fulfilled: $M_{x}=0$, $M_{y}=0$, or $M_{x}M_{y}^{*}$ is non-zero but real, $S_{3}$ is zero and the light is fully linearly polarized. Note that the specific choice of the $x,y,z$ directions does not change the conclusion obtained from Eq. 15, i.e. the polarization of the emitted lightHecht (2002). In the case of anisotropic systems such as MBP nanostructures, by choosing the $y$ axis along the dark transition direction and the $x$ axis along the bright transition direction (namely ${\bf\it M}=(1,0,0)$), from Eq. 15 it can be immediately deduced that the emitted light is linearly polarized.

In contrast, in-plane isotropic 2D systems with valley degeneracy like \ceMoS2 may have possible normalized dipole moments ${\bf\it M}$ that give other types of polarization (circular or elliptical) due to the mixed exciton states between $K$ and $K^{\prime}$ valleys. For example, when the exciton has ${\bf\it M}=1/\sqrt{2}(1,i,0)$, we obtain $S_{3}/S_{0}=1$, and the light is fully circularly polarized. Mixed exciton states can go through a fast decoherence process that leads to a significant loss of polarizationChen et al. (2018). However, such decoherence mechanism does not exist in MBP due to its anisotropicity which may be advantageous for quantum information applications that require long coherence time of quantum states. The above conclusions agree with the experimental photoluminescence measurements of the first exciton in MBPQiao et al. (2014); Wang et al. (2015); Li et al. (2017), which show a nearly perfect linear polarization.

In addition, for anisotropic systems like the MBP nanostructures we found that the light emission intensity ($I/(n({\bf\it Q}_{\textrm{ex}})\Gamma_{0}$)) always has a $\cos^{2}(\varphi_{L})$ angle dependence, where $\varphi_{L}$ is the azimuth angle of photon wavevector. This is the case regardless of the dimensionality and the polarization of the laser that excites the system. This behavior is different from in-plane isotropic 2D systems like BN and MoS2 that can have spherical symmetric or $\cos^{2}(\varphi_{L})$ angle dependence based on the mixture between $K$ and $K^{\prime}$ valleys (polar plots of photoluminescence for different 2D systems can be found in SI).

In conclusion, we presented a general framework to study exciton radiative lifetimes and light emission intensities in 2D, 1D, and 0D systems from many body perturbation theory. Based on it, important insights were provided on dimensionality and anisotropicity effects on exciton dynamics in phosporene nanostructures. For each dimensionality reduction the energy and momentum conservation leads to an additional factor of the order of $10^{-3}$ in the decay rates. This explained the general experimental observations that the exciton radiative lifetime is much longer in lower dimensional systems. Furthermore, in order to obtain the accurate radiative lifetime at finite temperature and quantitatively compare with experiments, accurate exciton energy dispersion in momentum space beyond the effective mass approximation must be considered. Finally, we demonstrated that the exciton anisotropicity in MBP nanostructures always leads to linearly polarized light emission unlike in-plane isotropic 2D materials such as \ceMoS2 and BN, which can emit light with arbitrary polarization.

Acknowledgement

We thank Francesco Sottile and Matteo Gatti for helpful discussions. This work is supported by NSF award DMR-1760260 and Hellman Fellowship. D.R. acknowledges financial support from Agence Nationale de la Recherche (France) under Grant No. ANR-15-CE29-0003-01. This research used resources of the Center for Functional Nanomaterials, which is a US DOE Office of Science Facility, and the Scientific Data and Computing center, a component of the Computational Science Initiative, at Brookhaven National Laboratory under Contract No. DE-SC0012704, the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the US Department of Energy under Contract No. DEAC02-05CH11231, the Extreme Science and Engineering Discovery Environment (XSEDE) through allocation TG-DMR160106, which is supported by National Science Foundation Grant No. ACI-1548562Towns et al. (2014).