Determining the directional strain shift coefficients for tensile Ge: a combined x-ray diffraction and Raman spectroscopy study

Until recently, strain engineering in the SiGe material system mainly concentrated on the improvement of electronic properties in tensile strained Si and compressively strained SiGe regions. Extending the approach to tensile strained Ge might enable a major breakthrough for the optical properties: tensile strain is predicted to revert the band structure of Ge from indirect to direct gap transitions which allows light emission and brings optical components integrated with Si technology within reach. Different concepts have been used to realize tensile Ge [1–3], often based on microstructuring. For fast and reliable device design and optimization, determining the strain distribution on a micron or even sub-micron scale is essential.

While x-ray diffraction (XRD) is the method of choice to quantitatively determine the strain state of crystalline matter, the experiments get elaborate for studies with sub-micrometer resolution, which in most cases require synchrotron radiation. An alternative, frequently used method to study strain locally and fast is micro-Raman spectroscopy (RS), which offers high spatial resolution, is non-invasive like XRD, available in many research labs and comparatively simple to set up [4, 5]. Raman spectroscopy, however, measures strain indirectly via the strain-induced shift of Raman-active phonon modes. This requires accurate knowledge of the strain shift coefficient b, which relates the Raman frequency shift $\Delta \omega$ to the absolute strain variation $\Delta \varepsilon$ inducing the shift [6]:

$$\begin{eqnarray}&&b=\frac{ \Delta \omega}{ \Delta \varepsilon}\end{eqnarray} \tag{ 1 }$$

b has been extensively studied for Si, Ge and SiGe alloys in the compressive strain region (e.g. [6–10]). However, in case of tensile Ge there is still no common ground mainly due to the fact that large tensile strains are not easy to realize technologically and an in situ calibration is technically challenging.

The variation of b with strain can be calculated from phonon deformation potentials (PDPs) p, q and r (see below and for example [11]). But literature values for phonon deformation potentials as well as b differ appreciably from each other, which makes strain measurements with RS arduous. For tensile strained Ge the reported PDPs are summarized in table 1.

Table 1. List of reported phonon deformation potentials in literature.

Authors	$p/{{\omega}_{0}}$	$q/{{\omega}_{0}}$	$r/{{\omega}_{0}}$	Year
Cerdeira et al [7]	−1.473	−1.933	$-{{1.087}^{\text{a}}}$	1972
Lockwood et al [12]	−1.45	−1.95	−1.1	1991
Peng et al [13]	−0.2	−1.25	−3.5	2009
Reparaz et al [6]	−1.66	−2.19	−1.11	2008

^aIn the original paper the value of  −10.87 is given. However, we believe that a decimal point fault could have occurred, since such a high value is not very likely; for comparison, the value for Si, which should be the same order of magnitude, is  −0.67 [7].

In the presented study the calibration of the Raman strain shift coefficients for tensile strained Ge is undertaken, by combining nano-focused high resolution XRD and ex situ power-dependent RS, together with a novel approach for fitting PDPs.

The strain shift coefficient defined in equation (1) is actually a directional quantity, i.e. it depends on the crystallographic orientation in which the strain is applied. To obtain a complete picture, shifts have to be measured for strains along different crystallographic orientations, which then allows to determine the phonon deformation potentials (once they are all known, strain shift coefficients for arbitrary strain states can be calculated).

While large multiaxial tensile strains are technologically extremely difficult to realize, large uniaxial strains have been achieved with different approaches [3, 14–16]. For this study, bridge structures as described below and in [3] with uniaxial strain along various crystallographic directions have been utilized. Scanning electron microscopy images are shown in figure 1. The investigation of these strains allows to extract the directional dependence of the strain shift coefficient.

Zoom In Zoom Out Reset image size

As described in [17] there are three Raman-active optical phonon modes for Si and Ge: two transverse optical and one longitudinal optical phonon mode. They are degenerate at k  =  0 due to the diamond cubic lattice structure.

If non-hydrostatic strain is applied to the crystal, this cubic symmetry is destroyed, which lifts the degeneracy. The splitting in phonon frequencies and the mixing of the phonon modes then contain complete information on the stress which is applied to the crystal. The total intensity I of the different phonon modes, for a given scattering geometry, can be determined by examining the Raman polarizability tensors ${{\mathbf{R}}_{j}}$:

$$\begin{eqnarray}&&I=C\underset{j}{\sum}\,|{{\mathbf{e}}_{i}}\centerdot {{\mathbf{R}}_{j}}\centerdot {{\mathbf{e}}_{s}}{{|}^{2}}\end{eqnarray} \tag{ 2 }$$

This is also called the polarization selection rule [11], where C is a constant related to the Raman scattering cross section, ${{\mathbf{e}}_{i}}$ the polarization of the incident light and ${{\mathbf{e}}_{s}}$ the polarization of the scattered light. The polarizability tensors ${{\mathbf{R}}_{j}}$ in coordinates along cubic $\langle 1\,0\,0\rangle$ directions with corresponding eigenvectors ${\bf \hat{e}}_j$ are given as:

$$\begin{eqnarray}{{\mathbf{R}}_{1}}=\left(\begin{array}{c} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{array}\right),~{{\mathbf{R}}_{2}}=\left(\begin{array}{c} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{array}\right),~{{\mathbf{R}}_{3}}=\left(\begin{array}{c} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{array}\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} {\bf \hat{e}}_1 = \left( \begin{array}{ccc} 1 & 0 & 0 \end{array} \right), ~ {\bf \hat{e}}_2 = \left( \begin{array}{ccc} 0 & 1 & 0 \end{array} \right), ~ {\bf \hat{e}}_3 =\left( \begin{array}{ccc} 0 & 0 & 1 \end{array} \right). \end{eqnarray} \tag{ 4 }$$

The effect of stress on the Raman modes for the three optical phonon modes are given by

$$\begin{eqnarray}&&\omega^2 {\bf \hat{e}}_j = \sum_k K_{jk}{\bf \hat{e}}_k\end{eqnarray} \tag{ 5 }$$

with j, k  =  1  −  3 [11, 18]. The ${\bf \hat{e}}_j$ are the eigenvectors in cubic coordinates, ω is the Raman mode frequency in the presence of strain and K_ij are the elements of the force constant tensor. These K_jk can be expanded in powers of the strain as follows:

$$\begin{eqnarray}&&{{K}_{jk}}=\omega _{0}^{2}{{\delta}_{jk}}+\underset{lm}{\sum}\,{{\varepsilon}_{lm}}K_{lmjk}^{\left(\varepsilon \right)}\end{eqnarray} \tag{ 6 }$$

${{\omega}_{0}}$ is the Raman frequency of the unstrained crystal, ${{\delta}_{jk}}$ the Kronecker delta and ${{\varepsilon}_{lm}}$ the elements of the strain tensor.

For Si and Ge, which crystallize both in diamond lattice structure, the symmetrical tensor $\mathbf{K}$ has only three independent and non-zero elements, which are also called the phonon deformation potentials:

$$\begin{eqnarray}\begin{array}{*{35}{l}} K_{1111}^{\left(\varepsilon \right)}=K_{2222}^{\left(\varepsilon \right)}=K_{3333}^{\left(\varepsilon \right)} & = & p \\ K_{1122}^{\left(\varepsilon \right)}=K_{1133}^{\left(\varepsilon \right)}=K_{2233}^{\left(\varepsilon \right)} & = & q \\ K_{1212}^{\left(\varepsilon \right)}=K_{1313}^{\left(\varepsilon \right)}=K_{2323}^{\left(\varepsilon \right)} & = & r \end{array}\end{eqnarray} \tag{ 7 }$$

From equations (5)–(7) the following secular equation follows:

$$\begin{eqnarray}\left|\begin{array}{c} p{{\varepsilon}_{11}}+q\left({{\varepsilon}_{22}}+{{\varepsilon}_{33}}\right)-\lambda & 2r{{\varepsilon}_{12}} & 2r{{\varepsilon}_{13}} \\ 2r{{\varepsilon}_{12}} & p{{\varepsilon}_{22}}+q\left({{\varepsilon}_{33}}+{{\varepsilon}_{11}}\right)-\lambda & 2r{{\varepsilon}_{23}} \\ 2r{{\varepsilon}_{13}} & 2r{{\varepsilon}_{23}} & p{{\varepsilon}_{33}}+q\left({{\varepsilon}_{11}}+{{\varepsilon}_{22}}\right)-\lambda \end{array}\right|=0\end{eqnarray} \tag{ 8 }$$

The eigenvalues ${{\lambda}_{k}}$ are related to the shifted phonon frequencies ${{\omega}_{k}}$ as

$$\begin{eqnarray}&&{{\lambda}_{k}}=\omega _{k}^{2}-\omega _{0}^{2}=\left({{\omega}_{k}}+{{\omega}_{0}}\right)\left({{\omega}_{k}}-{{\omega}_{0}}\right)\approx 2{{\omega}_{0}}\left({{\omega}_{k}}-{{\omega}_{0}}\right)\end{eqnarray} \tag{ 9 }$$

which gives

$$\begin{eqnarray} &&\Delta {{\omega}_{k}}={{\omega}_{k}}-{{\omega}_{0}}\approx \frac{{{\lambda}_{k}}}{2{{\omega}_{0}}}.\end{eqnarray} \tag{ 10 }$$

Solving the secular equation (8) gives then a relation between the Raman frequency of the three modes as a function of the strain [11].

This leads in general to quite complicated expressions, but as soon as the assumption of uniaxial stress holds, which is the case for the investigated structures, simplified expressions are determined. An elaborate finite element method (FEM) evaluation to proof this is shown in the appendix.

Thus, solving the equations above in the coordinate systems of the strained structure, i.e. the main axis in the direction of the uniaxial strain, simplifies the calculations (see [11]). A listing of the required steps rotating the different matrices and tensors into the proper orientations is given below.

The Raman shift and intensity variation induced by the strain was measured for different directions. In order to calibrate the strain shift coefficients and fit the PDP an absolute strain reference is needed. XRD is capable to determine the absolute strain value non-destructively and thus is the method of choice. Since the part of the bridges for which the uniaxial strain assumption is valid is only several μm large (see figure 1 and the appendix), it is important that both Raman and XRD probe only this part of strained material. This is ensured using a μ-Raman setup with a sub-micron spot size. In order to probe a similarly small volume also with the x-rays, synchrotron radiation with a nanofocusing optics was employed. Details of both setups are given below.

The strain components determined by XRD were then used together with the measured Raman strain shifts $\Delta \omega$ to fit the angle-dependent strain shift coefficients b in order to determine the PDPs for the tensile strain region.

Figure 1 shows scanning electron microscopy (SEM) images of the investigated structures. They are based on a Ge layer grown epitaxially on a SOI substrate with [0 0 1] surface orientation. The layer relaxes plastically at growth temperature, and during cooling to room temperature, tensile biaxial in-plane strain of 0.15% builds up due to different thermal expansion coefficients of Ge and Si. To achieve high strain values, dumbbell-shaped 'bridge' structures are etched out of the Ge layers, so that the strain is concentrated in the thin bridges, where it reaches values up to 0.8% in our samples⁵ (see below). To measure the strain shift coefficient in different crystallographic directions, bridges oriented differently with respect to the primitive crystal axes were fabricated, as indicated in figure 1. The geometry of the bridges gives actually rise to uniaxial stress along the constriction of the bridges, which was also confirmed in FEM simulations (see appendix). Details on the sample preparation can be found in [3] and [20].

The Raman strain shift and intensity has been measured on a WITec CRM200 confocal Raman system with a 532 nm excitation laser. Data has been recorded using a 100×  objective lens with a numerical aperture of 0.9 to couple light into a grating spectrometer (2400 lines mm⁻¹) via a multimode fiber with core and cladding diameters of 50 and 125 μm, respectively. The polarizer and analyzer setting was parallel to the bridge direction, corresponding to a Z(X,X)$\underline{\rm Z}$ geometry in Porto's notation. More experimental details on the setup can be found in [5]. With this setup it is possible to map the Raman shift over μm ranges. An optical microscope is used to align the excitation laser with a spot size of 720 nm to the centre of the bridge. The Raman shifts for differently oriented bridges were obtained with respect to the laser power and interpolated to zero laser power to yield the true Raman shifts, i.e. eliminate any influence of the local heating by the laser on the strain state. The resulting Raman shifts as a function of bridge orientation are shown in figure 2.

Zoom In Zoom Out Reset image size

A trend of decreasing Raman shifts with larger bridge inclination with respect to the [1 0 0] axis is observed, consistent with the variation of the strain shift coefficients with crystallographic orientation [7].

In order to determine the strain at the centre of the bridges directly, an x-ray nanodiffraction experiment has been conducted at beamline ID13 at the ESRF (details on the beamline can be found in [21]). The x-ray beam at an energy of 15.198 keV was focused with compound refractive lenses to $600\times 200$ nm² and the angular range around the nominal Ge Bragg peak was scanned while the individual bridges were moved through the beam employing a piezo-electric scanning stage mounted on the diffractometer The 2D maps in real space were obtained by scanning the sample position for every angular position [20].

At every real space position the actual (0 0 8) and (2 2 8) Bragg peak position is determined by calculating the centre of mass (COM). Possible tilt corrections obtained from the symmetric reflection are applied to the asymmetric reflection. From the COM in Q-space the in-plane and out of plane lattice constant and hence the in-plane and out of plane strain was determined.

The detailed strain variation along the bridge for bridge G-1 is shown in figure 3.

Zoom In Zoom Out Reset image size

The error bars are only shown for several datapoints for clarity. They contain the experimental uncertainty due to the angular goniometer resolution, the broadening of Bragg peaks due to a small tilt distribution of the bridges within the illuminated area, and due to the finite divergence of the focused beam. The latter effect actually dominates the experimental uncertainty, but has to be viewed as a worst-case estimate: the primary beam profile remains the same for all data points, and the COM determination of the peak position systematically determines the same center, hence the indicated error bars overestimate the true 'point-to-point' uncertainty.

To obtain the phonon deformation potentials from the experimental data, we employed the following scheme:

• The stiffness tensor $\mathbf{C}$ and strain matrix $\varepsilon$ are constructed in the coordinate system of the crystal, i.e. ${\bf \hat{e}}_i$ along cubic $\langle 1\,0\,0\rangle$ directions.
• In order to solve Hooke's law, i.e. retrieve strain components depending on the applied stress, both $\mathbf{C}$ and $\varepsilon$ are rotated into the coordinate system of the bridge. Therefore, tensor and matrix transformation rules are applied, which gives ${{\mathbf{C}}^{\prime}}$ and ${{\varepsilon}^{\prime}}$, respectively.
• In this bridge geometry uniaxial stress holds, which means the corresponding stress matrix ${{\sigma}^{\prime}}$ has only one non-zero element.
• Hooke's law is solved in this geometry, which is then a set of linear equations in seven variables (6 strain components and 1 stress component) and can be solved as a function of only one variable. The one unassigned variable in the resulting strain tensor ${{\varepsilon}^{\prime}}$ is chosen to be $\varepsilon _{11}^{\prime}$, which is the strain along the bridge.
• The strain tensor ${{\varepsilon}^{\prime}}$ is then rotated back into the crystal coordinate system and the components of this tensor $\varepsilon$ are used to solve the secular equation in crystal coordinates for eigenvalues ${{\lambda}_{k}}$ and eigenvectors $\mathbf{e}_{k}^{\text{strain}}$.
• The new Raman tensors $\mathbf{R}_{k}^{\text{strain}}$ in presence of strain are calculated as [17]
  $$\begin{eqnarray}&&\mathbf{R}_{k}^{\text{strain}}=\underset{j}{\sum}\,{{\mathbf{R}}_{j}}\mathbf{e}_{k}^{\text{strain}\{j\}},\end{eqnarray} \tag{ 11 }$$
  where $\mathbf{e}_{k}^{\text{strain}\{j\}}$ is the j^th component of the k^th eigenvector.
• The Raman intensities of the different modes are calculated according to [11]:
  $$\begin{eqnarray}&&{{I}_{k}}=\,|{{\mathbf{e}}_{i}}\centerdot \mathbf{R}_{k}^{\text{strain}}\centerdot {{\mathbf{e}}_{s}}{{|}^{2}}\end{eqnarray} \tag{ 12 }$$
• The average Raman shift is then:
  $$\begin{eqnarray} &&\Delta \omega =\frac{\underset{k}{\sum}\,{{I}_{k}}{{\lambda}_{k}}}{\underset{k}{\sum}\,{{I}_{k}}}\end{eqnarray} \tag{ 13 }$$
• This Raman shift enters into the strain shift coefficient given in equation (1), with $\Delta \varepsilon$ determined from the XRD measurements. The experimental XRD values along the bridge have been averaged over several data points along the centre of the bridge.
• A numerical fit of b is performed in order to determine the values for the deformation potentials, considering all different orientation angles β of the bridges with respect to the $\left[1\,0\,0\right]$ direction (see figure 1). The numerical fit generally depends on p, q, and r. However, r can not be determined in this case due to the uniaxial stress along the bridges. This implies that the resulting shear-strains ${{\varepsilon}_{ij}}$ with $i\ne j$ are virtually zero (discussed in the FEM section of the appendix) and r drops out in the secular equation (8). Therefore, it is only possible to fit p and q.

By combining the obtained Raman strain shifts and XRD strain values into the strain shift coefficient b, the optimum phonon deformation potentials have been fitted.

Using the literature values from [7] as starting point for the fit, the obtained values for $p/{{\omega}_{0}}$ and $q/{{\omega}_{0}}$ are $-1.566\pm 0.553$ and $-1.716\pm 0.175$, respectively, which gives the directional dependence of the strain shift coefficient b as plotted in figure 4.

Zoom In Zoom Out Reset image size

However, using the values from [6] or [13] as starting point, gives slightly different fitting results for p and q. The reason becomes clear from inspecting the fitting residuum over the parameter space closer. Figure 5 shows the residuum values over a reasonable range of the parameters p and q. It is found that the same minimum value is obtained for all pairs of p and q along a line

$$\begin{eqnarray}&&p=\left(2.86\pm 0.02\right)\centerdot q+\left(3.33\pm 0.02\right).\end{eqnarray} \tag{ 14 }$$

Zoom In Zoom Out Reset image size

Which particular minimum along this line is found by the fitting algorithm depends on the starting point used and the particular algorithm chosen.

Figure 5 also shows the literature values listed in table 1. Most of them show some deviations from our fit result, except for the values from [13], which fulfills quite exactly the linear relation in equation (14).

The consequence of using different deformation potentials is depicted in a different representation in figure 6. Here, the strain shift coefficients are plotted in a polar coordinate system as a function of the in-plane angle for different surface normal directions: The [0 0 1] case investigated in this study is depicted in red, but using the found deformation potentials, also other cases like [1 1 0] surface orientation (gray) and [1 1 1] surface orientation (blue) can be calculated, and are also shown in the plot. The lines represent the strain shift coefficients based on our fit result, where we have used the literature value from [7] for r. The shaded areas show the range of results obtained from the different data sets in literature. Our results show that present literature values rather overestimate the strain shift coefficients, leading to an underestimation of the actual strain when converting Raman line shifts.

Zoom In Zoom Out Reset image size

With this experiment it was possible to obtain the directional strain shift coefficient b for Ge in the tensile strain region. Experimentally determined Raman strain shifts and absolute XRD strain values have been used in a novel fitting approach to fit PDPs. Thus b could be determined for all directions ($\beta =0$–${{45}^{\circ}}$, other directions are symmetry-equivalent in the case of [0 0 1] surface orientation and uniaxial strain). The fitted phonon deformation potentials p and q are shown to be linearly dependent in the form of $p=2.858\centerdot q+3.333$. Deformation potential values for Ge from literature lead to overestimation of strain shift coefficients, and hence to underestimation of strain values obtained from Raman studies. So far, the available sample set did not allow us to determine all deformation potentials of Ge independent of each other. To do so and to improve the accuracy, the combination of nanobeam XRD and micro-Raman can be used with modified samples: Using different surface orientations will reduce the numerical fitting errors; in order to determine all three deformation potentials, samples with different strain states, especially also including finite shear strains, are required. The method can, however, easily be used to calibrate Raman-based strain measurements also in other material systems.

The work has been supported by FWF, Vienna (SFB IR-On, F2507-N08) and the Swiss National Science Foundation (SNF project no. 130181). The authors acknowledge ScopeM, ETH Zürich for the use of their facilities and support. Andreas Wyss acknowledges funding by the Helmholtz Gemeinschaft in the form of Helmholtz Virtual Institute VI530. The authors thank Manfred Burghammer, scientist in charge of beamline ID13 at the ESRF in Grenoble, for his help during the XRD experiment.

In order to verify the results obtained from Raman spectroscopy and XRD the bridges were also modeled in COMSOL Multiphysics, a finite element method (FEM) software. The bridges were all constructed with the dimensions obtained from SEM images as shown in figure 1. The model consisted of a 1.5 μm thick Ge layer with a tensile in-plane pre-strain of 0.15%. This pre-strain has been obtained from a reference sample with a conventional large-beam laboratory XRD measurement, being primarily sensitive to unprocessed areas between the bridges.

This pre-strained Ge layer follows a 1 μm thick layer of SiO₂ on top of the thick Si handle wafer. The used elastic constants are given in table A1. The FEM models were calculated under the assumption of linear elasticity. The boundary conditions were set to fixed constraints at the bottom of the structure and to roller conditions at the side walls framing the model. I.e. the structure was fixed at the bottom and could move only in vertical direction.

An exemplary model of the bridge structure A-1 is depicted in figure A1. The bridge is oriented at $\beta ={{27}^{\circ}}$ with respect to the [1 0 0] direction (x-axis). The contour-plot represents the calculated in-plane strain component along the bridge ${{\varepsilon}_{\parallel}}$.

Zoom In Zoom Out Reset image size

The procedure was repeated for all bridge orientations, and in each case the strain variations along each bridge were determined. Figure A2 shows the results of ${{\varepsilon}_{\parallel}}$ and ${{\varepsilon}_{\bot}}$ as well as the shear components ${{\varepsilon}_{xz}}$ and ${{\varepsilon}_{yz}}$ of all bridges.

Zoom In Zoom Out Reset image size

It can be seen that the in-plane strain ${{\varepsilon}_{\parallel}}$ decreases with increasing inclination with respect to the [1 1 0] direction. At the same time, also the absolute value of the out-of-plane strain ${{\varepsilon}_{\bot}}$ decreases with increasing inclination. To confirm the assumptions of virtually zero shear components, also ${{\varepsilon}_{xz}}$ and ${{\varepsilon}_{yz}}$ are plotted. The little spikes at the bridge ends and the position where the under-etching stops arise from the finite mesh size of the FEM simulation and actually overestimates the real values of the shear components. So the mentioned assumption of zero shear components holds.

Figure A3 shows the inplane and out-of-plane strains determined both from XRD and FEM at the centre of all the different bridges. It can be seen that the general trend in the experimental data follows very well the values obtained from FEM. The small quantitative differences between the absolute strain values in experiment and simulation most probably arise due to small differences between the idealized FEM model and the true sample geometries. For instance, the very critical length ratio of the wide underetched region and the thin constrictions might be slightly wrong due to non-perfect underetching of the Ge layer at the very rim of the pattern.

Zoom In Zoom Out Reset image size