Strain broadening of the 1042-nm zero-phonon line of the NV- center in diamond: a promising spectroscopic tool for defect tomography

TL;DR

This paper investigates strain broadening of the NV- center's IR zero-phonon line in diamond, developing a model to interpret variations with temperature and photoexcitation, and proposing its use for defect tomography.

Contribution

It introduces a model for strain broadening effects in the NV- center's IR ZPL, revealing new insights into defect contributions and temperature effects, and suggests a novel spectroscopic tool for defect analysis.

Findings

Identification of anomalous temperature-dependent broadening

Detection of defects other than nitrogen contributing to strain

Proposal of IR ZPL as a defect tomography tool

Abstract

The negatively charged nitrogen-vacancy (NV-) center in diamond is a promising candidate for many quantum applications. Here, we examine the splitting and broadening of the center's infrared (IR) zero-phonon line (ZPL). We develop a model for these effects that accounts for the strain induced by photo-dependent microscopic distributions of defects. We apply this model to interpret observed variations of the IR ZPL shape with temperature and photoexcitation conditions. We identify an anomalous temperature dependent broadening mechanism and that defects other than the substitutional nitrogen center significantly contribute to strain broadening. The former conclusion suggests the presence of a strong Jahn-Teller effect in the center's singlet levels and the latter indicates that major sources of broadening are yet to be identified. These conclusions have important implications for the understanding of the center and the engineering of diamond quantum devices. Finally, we propose that the IR ZPL can be used as a sensitive spectroscopic tool for probing microscopic strain fields and performing defect tomography.

Tables

Table 1. Table 1: The volume strain per unit concentration of various defects in diamond as we have calculated using DFT and as obtained previously lang91 by experiment.

$V\left[i\right]/V$	N0	N+	NV-	NV0
Calc.	0.23	0.01	$-0.16$	0.00
Expt.	0.41	–	–	–

Equations

$$\begin{array}[]{c}\displaystyle\nu_{\pm}=V(A_{1})e(A_{1})+V(A_{1}^{{}^{\prime}})e_{1}(A^{{}^{\prime}})\pm\\ \\ \displaystyle\sqrt{V^{2}(E)[e_{1}(E)+e_{2}(E)]^{2}+V^{2}(E^{{}^{\prime}})[e_{1}(E^{{}^{\prime}})+e_{2}(E^{{}^{\prime}})]^{2}}\>.\end{array}$$

$$g(\text{e})=\frac{3\gamma }{4{\pi }^2}[e_1(E{)}^2+e_2(E{)}^2+e_1(E^{{}^{\prime }}{)}^2+e_2(E^{{}^{\prime }}{)}^2+{\gamma }^2{]}^{-5/2},$$

$$\begin{array}[]{c}\displaystyle I(\nu)\propto\int d\textbf{e}g(e)(exp\left\{-\frac{[\nu-\nu_{0}-\nu_{-}]^{2}}{2\Delta^{2}}\right\}+\\ \\ \displaystyle exp\left\{-\frac{[\nu-\nu_{0}-\nu_{+}]^{2}}{2\Delta^{2}}\right\})\>,\end{array}$$

$$c[{\mathrm{N}}^0](P)=c[{\mathrm{N}}^0](0)\frac{1}{1+\Omega P}$$

$$\gamma (P)=\{\gamma [\mathrm{other}{]}^2+[\gamma [{\mathrm{N}}^0](0)\frac{1}{1+\Omega P}{]}^2{\}}^{1/2}$$

Full text

PACS numbers

61.72.Hh, 71.55.-i, 76.70.Hb, 61.72.jd

††preprint: APS/123-QED

I Introduction

In recent years, the negatively charged nitrogen-vacancy (NV-) center in diamond has attracted numerous research efforts aimed towards a wide range of applications. This color center exhibits long ground-state spin coherence and optically induced spin polarization over a broad range of wavelengths in ambient conditions, which makes it a perfect candidate for quantum sensing applications,qs1 ; qs2 ; qs3 ; qs4 ; qs5 ; qs6 room temperature operated masers,maser and quantum computing.qc1 ; qc2 ; qc3

The inhomogeneous broadening of the optical and spin resonances of NV- centers is a key limitation of many of these applications. It is currently believed that this inhomogeneous broadening is principally caused by the inhomogeneous strain and electromagnetic fields generated by random microscopic distributions of defects in the diamond, such as the substitutional nitrogen (N) center.stoneham69 ; davies70 ; davies71 However, the correlation between microscopic defect density and inhomogeneous broadening has not yet been quantitatively studied, and so this hypothesis has not yet been rigorously tested. Better understanding of the sources of inhomogeneous broadening will ultimately support superior engineering of diamond quantum devices.

The NV- center manifests two zero-phonon lines (ZPL) [Fig.1 (a)], the visible 1.946 eV (637 nm) and infrared (IR) 1.190 eV (1042 nm) associated with triplet-triplet ${{}^{3}}\textit{A}{{}_{2}}\leftrightarrow{{}^{3}}\textit{E}$ and singlet-singlet ${{}^{1}}\textit{E}\leftrightarrow{{}^{1}}\textit{A}{{}_{1}}$ transitions, respectively.Doherty_PhysRep_2013 The ${{}^{1}}\textit{E}\leftrightarrow{{}^{1}}\textit{A}{{}_{1}}$ transition plays a crucial role in the spin dynamics of the NV- center because, along with the triplet-singlet non-radiative intersystem crossings (ISCs) [shown by the dashed arrows in Fig. 1(a)], it is part of the spin-dependent channel that enables the optical spin polarization and readout of the center’s ground state spin.Doherty_PhysRep_2013 Thus, further study of the IR ZPL is important to the understanding of the properties of the singlet levels and the ISCs.goldman1 ; goldman2

In this paper we report a detailed study of the splitting and inhomogeneous broadening of the NV- center’s IR ZPL. We interestingly find that the splitting and broadening can be reduced by increasing excitation intensity or reducing excitation wavelength. To explain these observations, we develop a model for the inhomogeneous broadening that includes the strain induced by photo-dependent microscopic distributions of point defects. We use this model to distinguish strain broadening from other sources of broadening and, by doing so, identify an anomalous temperature dependent broadening mechanism, which potentially provides additional evidence of a strong Jahn-Teller effect in the 1E level. We show that the quantified strain broadening cannot be completely accounted for by the presence of N centers, thereby suggesting that there are major sources of broadening yet to be identified. We propose that the IR ZPL has the potential to be used as a sensitive spectroscopic tool to probe the microscopic distribution of strain fields and, with further work that identifies the major sources of strain, to perform tomography of defect concentrations.

II Experimental Method

Two different diamond samples containing NV- centers were studied in this work. The main difference between the samples was the concentration of N centers. Sample 1 (S1) was a type IIa single crystal diamond plate with dimensions 3$\times$3$\times$0.3 mm that was fabricated commercially by Element-6 using chemical-vapor deposition (CVD). The initial concentration of N was $\lesssim$100 ppb. The sample was irradiated to a total dose of 1018 cm*-2* with 900 keV electrons at room temperature. Sample 2 (S2) was a type Ib single crystal diamond plate with dimensions 3$\times$3$\times$0.3 mm that was fabricated commercially by Element-6 using high-pressure high-temperature (HPHT) growth and initially contained $\lesssim$200 ppm of N. The sample was irradiated with 3 - 5 MeV neutrons at temperatures close to room temperature to a total dose of 1018 cm*-2*. S2 was then annealed at a temperature T = 800∘C for two hours in the presence of a forming gas (H2) in order to achieve a high density of NV- centers via the capture of mobile vacancies by N centers.Mita_1996 ; Budker_PRB_NV_production_2010 S1 was not subjected to any thermal treatment. Thus, in summary, S1 has low densities of NV and N, whereas S2 has high densities.

The optical properties of the NV- centers in each sample was studied by performing micro-photoluminescence (PL) measurements using a Horiba Jobin-Yvon T64000 spectrometer equipped with a confocal microscope and a closed-cycle helium system for cryogenic microscopy. A Nd:YAG-laser line 2.33 eV (532 nm) and a He-Ne-laser line 1.96 eV (632.8 nm) were used for photoexcitation. The beam imaging, normal to the surface, was focused by a Mitutoyo 100xNIR (NA = 0.50) objective lens into each sample with a spot size of $\sim$2 $\mu$m. The same objective collected the PL. We used a CCD camera together with 1800 lines/mm and 600 lines/mm gratings to measure the PL spectra with the spectral resolution of $\approx$0.05 meV in the visible and IR ranges, respectively.

III Results and Analysis

We start with the comparison of the PL spectra of the NV- centers in samples S1 and S2 that are presented in Fig. 1(b,c). At a glance, there are significant differences in both the visible and IR ZPLs of the samples. In S1, narrow ZPLs were observed with simple Gaussian inhomogeneous lineshapes characterized by full widths at half maximum (FWHM) of $\sim$0.14 meV (visible) and $\sim$0.2 meV (IR). In S2, the visible and IR ZPLs differ significantly. The visible ZPL in S2 has shifted and broadened relative to S1 [inset in Fig 1(b)], whereas the IR ZPL in S2 has instead split and broadened.

The differing inhomogeneous lineshapes of the visible and IR ZPLs in S2 can be qualitatively explained by recalling that the ZPLs have different susceptibilitiesRogers_Doherty_2015 ; Doherty_NJP_2011 to inhomogeneous strain and electric fields induced by random distributions of defects.stoneham69 ; davies70 ; davies71 Being a type Ib diamond, the most prevalent defect in S2 is expected to be the N center. In thermal equilibrium, the N center adopts a neutral charge state N0. In this state, the center is paramagnetic and spontaneously undergoes a large distortion where one of the N-C bonds is substantially elongated.Loubser_vanWyk_1978 The known stress susceptibility parameters of the visible and IR ZPLsRogers_Doherty_2015 imply that the strain induced by a N0 center’s distortion will primarily shift the visible ZPL of a proximal NV- center, whilst primarily splitting its IR ZPL. Averaging over the strain fields induced by many N0 centers as well as over many NV- sites, one therefore expects to observe an inhomogeneously broadened visible ZPL without a splitting and an inhomogeneously broadened IR ZPL with a discernable splitting, as exhibited in Fig. 1(b,c).

Now, under photoexcitation with photon energy $\gtrsim 1.7$ eV,Farrer_1969 ; Isberg_PRB_2006 the N0 center can be photoionized to form the positive charge state N+. In this state, the center is non-paramagnetic and returns to an undistorted configuration where the N atom resides at the substitutional site.stoneham92 As a result, N+ generates an electric field due to its charge, but induces minimal strain. Since first-principles theory predicts that the IR ZPL is not susceptible to electric fields,Doherty_NJP_2011 the photoionization of N0 to N+ should lead to narrowing of the IR ZPL, thereby providing a means to test that N0 is indeed the dominant source of inhomogeneous broadening. The same is not true for the visible ZPL because it is known to be susceptible to electric fields.Doherty_PhysRep_2013 Therefore, any narrowing of the visible ZPL due to the smaller strain generated by N+ will be countered by additional broadening due to the electric field it produces. Consequently, in the following we focus on the IR ZPL as a tool to probe the sources of inhomogeneous broadening.

To quantitatively describe the inhomogeneous broadening of the IR ZPL, we adapt the model of defect-induced strain broadening proposed in Refs Klimin_Malkin_2010, and Malkin_PRB_2012, for rare earth doped inorganic crystals. The key assumptions of model are that the strain contributions from individual point defects add linearly and the crystal is well approximated by an isotropic elastic continuum. In the model, the strain-induced changes of the two transition energies $\nu_{\pm}$ of the IR ZPL are

[TABLE]

where e(A1) = exx + eyy + ezz, e(A1${}^{{}^{\prime}}$) = $\frac{1}{2}$(2ezz - exx - eyy), e1(E) = eyy - exx, e2(E) = 2exy, e1(E${}^{{}^{\prime}}$) = 2exz and e2(E${}^{{}^{\prime}}$) = 2eyz are the deformation tensor components that have been defined to have explicit C3v symmetry properties in the coordinate frame of the NV center, and V(A1), V(A1${}^{{}^{\prime}}$), V(E), V(E${}^{{}^{\prime}}$) are the corresponding strain susceptibility parameters of the IR ZPL. The susceptibility parameters can be obtained from those reported in terms of stress in the crystal coordinate system in Ref Rogers_Doherty_2015, by applying the elastic constants of the crystal and performing the required coordinate frame rotations. One finds that V(A1)= 636.48, V(A1${}^{{}^{\prime}}$)= –437.76, V(E)= –1567.17, V(E${}^{{}^{\prime}}$)= –1373.22 in units of meV/strain for the IR ZPL. Note that from the above expression for $\nu_{\pm}$, it can be seen that the IR ZPL splits due to the deformations of E symmetry, whilst it shifts due to those of A1 symmetry.

Following Ref. Malkin_PRB_2012, , assuming locally uniform densities of point defects, the statistical distribution function of the deformation tensor components e $\equiv$ [e1(E), e2(E), e${}_{1}($E${}^{{}^{\prime}}$), e2(E${}^{{}^{\prime}}$)] is

[TABLE]

where the distribution width $\gamma$ = \frac{\pi}{27}$$\frac{1+\sigma}{1-\sigma}$$\{\sum_{i}[c[i]\DeltaV$[i]/V)^{2}\}^{1/2}$ depends on the Poisson ratio of diamond $\sigma$ = 0.2,field12 the fractional concentration $c[i]$ of the $i^{th}$ species of point defect and the volume strain $\Delta V[i]/V$ per unit fractional concentration of defects of the $i^{th}$ species.Klimin_Malkin_2010 ; Malkin_PRB_2012 The strain broadened lineshape of the IR ZPL is then the convolution of g(e) and a Gaussian lineshape

[TABLE]

where the Gaussian width $\Delta$ accounts for the broadening induced by the deformations of A1 symmetry, all other inhomogeneous broadening (not arising from strain induced by point defects) and homogeneous broadening. In the above, the integration is taken over the four E symmetry deformation tensor components, and $\nu_{0}$ is the energy of the unperturbed IR ZPL. This integration must be performed numerically. Note that we have not explicitly treated the broadening induced by the deformations of A1 symmetry because, as supported by the following observations, the smaller strain shift of the IR ZPL implies that their contributions to the width $\Delta$ are negligible compared to the other sources of inhomogeneous and homogeneous broadening.

As demonstrated in Fig. 2(a), our expression for $I(\nu)$ successfully describes the inhomogeneous lineshape of the IR ZPL in S2. Fig. 2(b) depicts the temperature dependence of the IR ZPL under the same photoexcitation conditions. The dashed lines in Fig. 2(b) indicate that the ZPL splitting does not change with temperature, thereby implying that the strain induced broadening is temperature independent. The temperature dependence of the ZPL rather arises from the $\propto T^{3}$ temperature variation of the Gaussian width $\Delta$ that is plotted in Fig. 2(c). Such temperature dependent broadening has been observed in the NV- center’s visible ZPLFu_PRL_2009 and the ZPLs of other defects in diamond,jahnke15 and attributed to temperature dependent electron-phonon scattering. This is likely to be the case also for the IR ZPL, however it must be noted that some caution must be exercised because the Gaussian width $\Delta$ is not purely related to the homogeneous broadening that would arise from electron-phonon scattering, but also has contributions from other mechanisms of inhomogeneous broadening. If this complication is ignored, then our results suggest that the electron-phonon scattering affecting the IR ZPL has a somewhat unusual $\propto T^{3}$ temperature variation. Although a similar temperature variation has been observed for the ZPL of the negatively charged silicon-vacancy (SiV-) center, this has been shown to be a modification of the usual $\propto T^{5}$ Raman-type scattering mechanism by the large spin-orbit interaction of the SiV- center.jahnke15 No such spin-orbit interaction can be present in the singlet levels involved in the NV- center’s IR ZPL. It is possible that the $\propto T^{3}$ temperature variation of the IR ZPL arises from the presence of a strong Jahn-Teller effect,Hizhnyakov04 which is consistent with the dramatic asymmetry of photoabsorption and PL spectra of the IR ZPL’s vibrational sideband that is currently poorly understood.kehayias13 However, substantial further investigation is required before this conclusion can be drawn.

To test the role of N0 in strain broadening, we analysed the dependence of the IR ZPL on the photoexcitation energy and power. As can be seen in Fig. 3(a), increasing the photoexcitation energy (decreasing wavelength) with fixed power and temperature results in the decrease of $\gamma$ with almost constant $\Delta$. A result that is consistent over the two different positions in S2 that we sampled. Fig. 3(b) demonstrates that increasing the photoexcitation power with fixed energy and temperature similarly results in the decrease of $\gamma$ with almost constant $\Delta$. These results support the hypothesis described above that increased photoionization of N0 at higher photoexcitation energies and powers will reduce strain broadening. This hypothesis is further supported by the fit of the photoexcitation power dependence of $\gamma$ in Fig. 3(c) by a simple two-state rate equation model of N0 photoionization. Defining a power dependent photoionization rate $BP$ and a power independent recombination rate $A$, the steady state concentration of N0 is

[TABLE]

where $c[\mathrm{N}^{0}](0)$ is the concentration of N0 at zero power and $\Omega=B/A$. The steady state expression for $\gamma$ is then

[TABLE]

where $\gamma[\mathrm{N}^{0}](0)=\frac{\pi}{27}\frac{1+\sigma}{1-\sigma}c[\mathrm{N}^{0}](0)\frac{\Delta V[\mathrm{N}^{0}]}{V}$ is the strain broadening due to N0 at zero power and $\gamma[\mathrm{other}]$ is the power independent strain broadening due to all other sources. The parameter values obtained from the fit in Fig. 3(c) are $\gamma[\mathrm{other}]=1.30\pm 0.04$, $\gamma[\mathrm{N}^{0}](0)=1.27\pm 0.07$ and $\Omega=0.004\pm 0.002$ $\mathrm{\mu W}^{-1}$.

The power dependence of the strain broadening $\gamma(P)$ appears to imply the unexpected conclusion that N0 causes at most $\gamma[\mathrm{N}^{0}](0)/\{\gamma[\mathrm{other}]+\gamma[\mathrm{N}^{0}](0)\}\approx 50\%$ of the strain broadening, thereby indicating that major sources of strain broadening are yet to be identified. The other defects known to be in the sample are NV-, its neutral charge state NV0, and the vacancy (V) center. However, given that the concentrations of NV- and NV0 must be less than N0 by growth arguments, and that our annealing of S2 is expected to almost completely remove V centers, for these centers to be the major sources of strain broadening, they must induce much larger strains than N0.

Another issue arises in the quantitative analysis of the fit parameters. Using the value of $\Delta V[\mathrm{N}^{0}]/V=0.41$ that has been obtained from x-ray crystallography,lang91 the fitted value of $\gamma[\mathrm{N}^{0}](0)$ yields a N0 concentration of $c[\mathrm{N}^{0}](0)\approx 1775$ ppm at zero photoexcitation power, which is $\sim 10$ times larger than the specification of $\lesssim 200$ ppm that was assigned by the sample supplier. Possible origins of this overestimation are incorrect specification by the supplier, large non-uniformity of N0 concentration over the sample or photoionization of NV-. Considering the latter further, it is known that NV- photoionizes to NV0 under excitation with energy $>1.946$ eV.Doherty_PhysRep_2013 Thus, photoionization of NV- could occur under the excitation conditions we employed. However, similar to before, for photoionization of NV- to explain the overestimation, the difference in strain induced by NV- and NV0 must be much larger than the difference between N0 and N+.

To test these explanations, we have performed density functional theory (DFT) calculations of the strain induced by N0, N+, NV- and NV0 centers. Calculations were performed using the VASP programVASP and employed the screened hybrid functional of Heyd, Scuseria, and Ernzerhof HSE06 to describe the electronic structure. Interaction between ions and valence electrons were treated with the PAW approach.VASP Wavefunctions were expanded in plane waves (using 400 eV for the kinetic energy cutoff). We have used 216 (512) atom supercells to model N0 and N+ (NV- and NV0). The Brillouin zone was sampled at the $\Gamma$ point to make sure local symmetry is properly represented. Defect-induced stress was calculated as the average of the three diagonal stress components $p=(\sigma_{xx}+\sigma_{yy}+\sigma_{yy})/3$ in the supercell, mimicking the averaging over different defect orientations. Subsequently, $p$ was corrected for a small residual stress of the pristine supercell when calculated using exactly the same computational parameters (k-point sampling, cutoff energy, and the supercell size). Crucially, stress of charged defects (N+ and NV-) has been corrected for the spurious potential alignment term contribution using the procedure of Bruneval et al [Eq. (26) in the Ref. Bruneval_PRB_2015, ]. Absolute deformation potentials of bulk diamond states needed for that procedure have been taken from calculations reported in Ref. Li_PRB_2006, . Finally, volume expansion coefficient has been determined via $V[i]/V=p/(Bx)$, where $B=481$ GPa is the calculated bulk modulus of diamond (cf. the experimental value of 443 GPa), and $x$ is the concentration of defects in our supercells (either $1/216$ or $1/512$).

The results of our DFT calculations are summarized in Table 1. The value for the volume strain $(V[\mathrm{N}^{0}]/V)$ per unit fractional concentration of N0 underestimates the experimental value by \sim$$50\%. Part of this disagreement stems from the shortcomings of DFT and, likely, relatively large defect concentrations in actual calculations. However, the disagreement can also arise from experimental errors in determining the actual nitrogen concentration. In agreement with expectations, the value for N+ is negligible. The value for NV- is negative, which is consistent with the expectations based upon the presence of a vacancy in the defect. Interestingly, the value of $V[i]/V$ for NV0 is almost zero, which reflects a small outward relaxation of atoms adjacent to the vacancy in NV0 as compared to NV-. Importantly, the magnitudes of the values for N0 and NV- are comparable and thus, in accordance with the arguments above, we can conclude that NV- is not the unidentified major source of broadening nor is its photoionization the explanation for the overestimation of the N0 concentration. The resolution of these issues requires further investigation beyond this work. Such future work should include the confirmation of the average concentration of N0 in the samples using electron paramagnetic resonance and infrared absorption spectroscopy as well as the systematic study of large non-uniformities of the N0 concentration in the samples. Indeed, evidence of the latter can be immediately seen in the variation of $\gamma$ between the two positions in S2 that are depicted in Fig 3. If the issues can be resolved, then analysis of the strain broadening of the NV- center’s IR ZPL will be a promising tool for defect micro-/nanoscale tomography.

IV Conclusion

To summarize, we have reported a detailed study of the splitting and inhomogeneous broadening of the NV- center’s IR ZPL. We have identified that its inhomogeneous lineshape differs from the center’s visible ZPL because the IR ZPL has a different susceptibility to strain fields and is not susceptible to electric fields. This singular susceptibility motivated the use of the IR ZPL to quantify the broadening due to different microscopic distributions of defects via the strain fields that they generate. We developed a model of the inhomogeneous lineshape of the IR ZPL by adapting the strain broadening theory proposed in Refs Klimin_Malkin_2010, and Malkin_PRB_2012, . We applied the model to interpret observations of temperature variation of the IR ZPL shape. By distinguishing strain broadening from other sources of broadening, we identified an anomalous $\propto T^{3}$ temperature dependent broadening mechanism, which potentially provides additional evidence of a strong Jahn-Teller effect in the 1E level. This additional evidence furthers the understanding the singlet levels of the NV- center and their role in the center’s optical spin polarization and readout mechanisms.

We interestingly found that the strain broadening of the IR ZPL can be reduced by increasing excitation intensity or reducing excitation wavelength. Through the application of our model, these effects were found to be qualitatively consistent with the photodependence of the concentration of N0 centers. However, quantitative analysis revealed that only $\approx 50\%$ of the strain broadening can be attributed to the photodependence of N0, thereby indicating that there exist other, yet to be identified, sources of strain broadening. Furthermore, the quantitative analysis yielded an estimate of the local N0 concentration that was $\sim 10$ times larger than what the supplier specified as the average concentration. To investigate these issues further, we performed DFT calculations of the strain induced by different charge states of N and NV centers. These calculations supported our initial interpretation and thus did not resolve these issues of the quantitative analysis. We conclude that further work is required to identify the unidentified sources of broadening and to test the estimation of the local concentration of N0. Given this further work, we proposed that the IR ZPL has the potential to be used as a sensitive spectroscopic tool to probe the microscopic distribution of strain fields and to perform tomography of defect concentrations.

Acknowledgements.

This work has been supported by the ARC grants FT110100225, CE140100003 and DE170102735, and performed according to the Russian Government Program of Competitive Growth of Kazan Federal University. B. C. G acknowledges the support of the ARC Centre of Excellence for Nanoscale BioPhotonics (CNBP) and an ARC Future Fellowship. V. A. S. acknowledges the support of the RFBR Project 16-52-76017 (ERA.Net RUS Plus project ”DIABASE”). A. A. acknowledges the support of the Research Council of Lithuania via Grant No. M-ERA.NET-1/2015. Computations were performed at the High Performance Computing Center “HPC Sauletekis” in Faculty of Physics, Vilnius University. We wish to acknowledge fruitful discussions with Boris Malkin, Sergey Feofilov and Neil Manson.