Supplementary Information Supplementary Figures

Supplementary Figure 1. Optical and hyperspectral Raman images of an exfoliated black phosphorus flake with armchair and zigzag edges. (a) Optical microscope image of the measured flake; the white dashed rectangle defines the analysed area. The corresponding crystal orientation and crystallographic axes are presented in the inset. (b) Hyperspectral Raman intensity images of all observed modes (columns), for the XZ and ZX configurations (indicated on the left). Mode intensities are represented by the colorbar and were normalized to the highest value of each mode. Whenever necessary, for better visualization, the intensity was multiplied by a factor, as indicated at the bottom-right corner of each image. (c) Intensity profiles along the yellow lines shown in (a) for the A 1 g , A 2 g and B 2g modes. Columns indicate the different analysed modes and rows correspond to the parallel and orthogonal incident/analysed polarisation configurations. Orange and grey curves indicate z and x incident polarisation, respectively. Solid and dashed curves are the profiles taken at the zigzag and armchair edges, respectively.

Initial Relaxed  Figure 2 we compare, as an example, the ZZ spectra taken in one of the bright spots in Figure 1(d) of the manuscript and in a perfect region within the interior of the sample. In this configuration only the A 1 g and A 2 g are expected to appear. We can observe that the only difference in these two spectra is the presence on an intense background in the bright spot spectrum. As mentioned in the main manuscript text, this background is due to a luminescence process, showing that the bright spots might be related to sample imperfections and/or oxidation [1]. The important point is that the anomalous modes at the edges of the sample are not observed in the bright spots, which indicates that there is no relationship between the spectra at the edges and at the bright spots.

Results obtained in thinner samples
Measurements similar to those presented in Figure 1 of the main manuscript were made in two thinner flakes, with thicknesses of 30 nm and 6 nm. Supplementary Figure 3

(a)
shows hyperspectral Raman images for the A 1 g , A 2 g and B 2g modes in the ZZ, XX, XZ and ZX configuration, for the few-layer, 6 nm, flake. In the figure, the zigzag and armchair directions are along the horizontal and vertical axes, respectively. Despite the fact that the Raman signals are much weaker and that the Raman maps are much noisier than those of Fig.   1(d), we can observe the same anomalous results at the edges of the sample. In particular, we note the presence of the B 2g mode in the armchair edge for the ZZ configuration, while the A g modes appear in the armchair edge for the XZ configuration and in the zigzag edge for the ZX configuration. Also, Supplementary Figure 3(b) shows a comparison between the spectra of the three measured flake thicknesses, for an armchair edge in the XZ configuration.
Each spectrum was normalized by its B 2g mode intensity. The observation of the A g modes in the figure is indicative of the change in the atomic structure at the edges, which therefore continues to occur down to few-layer thicknesses. Notice that the amplitude of these modes decreases with crystal thickness. Besides, the non-expected modes, B 1g and B 1 3g , were not observed for the two thinner flakes since, even for the 300-nm thick flake, their intensities are very low.
In order to probe the temperature dependence of the Raman spectra at the edges and away from the edges (at the positions marked in Supplementary Figure 4(a)), the incident laser power was varied. All phonon frequencies were observed to downshift with increasing laser power. This is an expected behavior for materials with a positive thermal expansion coefficient. Supplementary Figure 4 Interestingly, the softening of the A 2 g phonon frequency with increasing laser power (temperature) is larger at the edges when compared to the interior of the sample, and the curve is steeper for the zigzag edge than for the armchair edge. We interpret these results as a consequence of the anisotropic thermal conductivity in BP [2]. Also, the differences in heat dissipation in the bulk and at edges lead to higher temperatures at the edges than in the sample interior.

Thickness dependence of the relative intensities of the BP and Si Raman peaks
In all measurements, the BP sample was on the top of a Si/SiO 2 substrate, allowing us to observe both the Raman signal of BP and of the Si substrate. Supplementary Figure 11 presents the spectra from flakes with 16 nm (black), 65 nm (blue) and 117 nm (green) thicknesses. We can clearly observe that the Si/BP relative intensities decrease with increasing thickness of the BP sample.

Crystallographic axis determination
It has been recently shown that the determination of the crystallographic axes of BP from the angular dependence of polarized Raman spectroscopy is not necessarily correct, since the Raman tensor elements depend both on the excitation laser energy and on flake thickness [3]. Therefore, the maxima and minima in the angular dependence can be associated with both the armchair or zigzag direction depending on the sample thickness. However, distinction between these two directions can be unambiguously made by the anisotropic optical absorption response of BP (linear dichroism) [4].
In this work, while obtaining the polarized Raman spectra of BP deposited on the top of a silicon substrate (with an oxide top layer), we also analysed the angular intensity dependence of the Si Raman peak at 525 cm −1 . First, with an analyser parallel to the polarization of the incident light, the BP crystal was rotated until the amplitude of its B 2g Raman band was null. At this orientation the polarization is aligned with either the armchair or zigzag direction [3]. Next, spectra were taken at the flake and away from it (just the substrate). The sample was then rotated by 90 • and another pair of spectra was taken. With the excitation wavelength at the visible range, the BP crystal is more absorbing along the armchair direction [5]. Therefore, the ratio between the silicon signal away from and through the flake is higher with polarization parallel to this direction. The solid lines in Supplementary Figure 9 show the Raman spectra through BP, while the dotted curves show the silicon spectrum measured away from the flake. Blue and red curves indicate the two principal orientations of the measured BP flake with a thickness of 170 nm. All spectra were taken using a 488-nm laser line and a 1-mW incident power. It is clear from the figure that, through BP, the silicon mode exhibits lower amplitude in the red curve, which is, thus, assigned to an armchair polarization. The reliability of this method was confirmed through comparison with electron diffraction and HRTEM data, which are presented in the next  Supplementary Figure 10(b) and (c), respectively. Notice that, to within experimental resolution, the electron diffraction patterns are the same in these two regions. The lattice parameters evaluated from the diffraction are a = 0.324nm and b = 0.424nm and are in agreement with previous results [6]. The diffraction patterns corroborate the edge determination obtained via our axis determination method.

Supplementary Note 2 -Theory and Methodology
For the theoretical calculations, plane-wave density functional theory was employed to obtain the electronic ground-state using the Perdew-Burke-Ernzerhof (PBE) [7] generalized gradient approximation exchange-correlation functional, currently implemented in the Quantum-Espresso package [8]. Van der Waals corrections within the semi-empirical dispersion scheme (PBE-D) proposed by Grimme [9] were included.
Norm-conserving pseudopotentials with 3s3p states were adopted to describe electronic states of phosphorus. The Brillouin zone was mapped with a k-sampling grid within the Monkhorst-Pack scheme using 7 × 7 × 7 for the bulk and a grid of 10 × 7 × 1 (1 × 8 × 8) for the armchair (zigzag) surfaces samples. The kinetic energy cutoff was set at 90 Ry and 100 Ry for bulk and slab geometries, respectively. Furthermore a vacuum region of 16Å was adopted for the supercell related to the surfaces. The structures were fully optimized to their equilibrium position with forces smaller than 0.002 eVÅ −1 and the unit cells were relaxed to a target pressure of 0.2 kbar. Supplementary Figure 5 shows the structures obtained after relaxation highlighting the unit cell used in the calculations. Supplementary Table 1 presents the unit cell lattice parameters before and after relaxation for bulk and for both edge terminations.
Within the linear response framework, the lattice distortion induced by phonons can be considered as a static perturbation acting upon the electrons. This enables the construction of the matrix elements for harmonic interatomic force constants where E is the electronic ground-state energy and u αi (R) represent the ith atom displacement in the αth direction at position R of the unit cell. By performing the Fourier transform of equation 1 one is able to solve the secular equation to obtain the phonon frequencies ω and normal modes w αi . Here the dynamical matrix D αi,βj (q) stands for the Fourier transform of the interatomic force constants.
The linear response approach can also be used to calculate the dielectric tensor (response which relates the bare electric field E 0 with the screened one E. In equation 3, E el is the electronic energy of the system in the presence of a uniform electric field E l(m) along direction l(m), and Ω is the system volume. As we will see bellow, the dielectric tensor is a crucial observable in determining the Raman intensities.
Within the Placzek approximation [10], the Raman intensity of a particular normal mode ν is calculated as: where e i (e s ) is the incident (scattered) light polarization, n ν is the Bose-Einstein distribution, and A ν is the Raman tensor with matrix elements: In equation 5, M γ corresponds to the atomic mass of the γ − th atom, E el is the electronic energy of the system, E l is the component of the electric field along direction l, u kγ is the displacement of the γth atom in the kth direction and w ν kγ is the orthonormal vibrational eigenmode ν [10]. Further details of the methodology are given in Ref [11].
The Raman tensor for all the structures were obtained using the PHonon code, currently integrated into the Quantum-Espresso package. It relies on the perturbative treating of electron-phonon interaction, in what is known as density functional perturbation theory (DFPT) [12].
For a quantitative analysis of the Raman intensities dependency on polarisation, we defined the polarisation vector for the incident light as e i = (sinθ 0 cosθ) T , the vector for scattered light in the same polarisation as e s = (sinθ 0 cosθ) T and the vector for scattered light at the orthogonal polarisation as e ⊥ s = (−cosθ 0 sinθ) T (the angle θ is measured with respect to z crystalographic axis, i.e., the zigzag direction).

Considering a generalized Raman tensor
one can obtain the following expressions for the measurable intensities as a function of the angle θ:

Bulk black phosphorus
In Supplementary Figure 6, we show the crystalline structure of AB-stacked black phosphorus as well as the atomic displacements of the six active Raman modes obtained within the DFPT formalism. The calculated Raman frequencies for the A 1 g , B 2g , and A 2 g modes are 343, 425 and 452 cm −1 , respectively. We also performed calculations using the AA stacking.
Our simulations, as well as others' [13], show that the AB stacking is a few meV per atom more stable. Furthermore, Supplementary Figure 13 shows the Raman spectra simulated using DFPT for both stackings compared to our experiments. There we can see that the ABstacked structure yields the best agreement with the experimental results (dashed vertical lines).
The angular dependence of the Raman spectra can be seen in Supplementary Figures 12(a) and 12(b) for the parallel and crossed polarisation configurations, respectively; it is evident that only modes A 1 g , B 2g , and A 2 g are observable, as predicted by group theory, because incidence is assumed to be along the y axis, as in the experiments.

Edge structures
For studying the edges of black phosphorus, we first calculated the atomic vibrational modes for both zigzag and armchair slabs by solving equation 2, and compared their intensities with those for bulk. We then looked for edge enhancements similar to those that were experimentally observed. Even though simulation-experiment agreement was not achieved for all modes in all scattering configurations, several experimental observations were qualitatively reproduced. Here, we highlight the situations for which such an agreement is achieved.
We begin our analysis by considering the zigzag case. There are 72 vibrational modes, which include confined, inactive, and active modes. We identify the Raman modes observed in the experiment by adopting two criteria: (i) we focus on a range of frequencies close to those calculated for the Raman active modes of the bulk material and (ii) we compare the atomic displacement vectors of the slab with those found for the bulk case. The results were presented and discussed in Figure 3 of the main manuscript.
In Supplementary Figure 7, we compare the Raman intensities of the bulk with those obtained simulating the zigzag slab, for the ZZ, XX, and XZ scattering configurations. The totally symmetric Raman modes for the slab arise at A 1 g = 344.6 cm −1 and A 2 g = 441.7 cm −1 , whereas the non-totally symmetric modes B 1g , B 1 3g , B 2g emerge at 175.4, 224.6 and 417 cm −1 respectively. Notice that the frequency related to the Aby 10 cm −1 with respect to the bulk, which might be associated to the quantum confinement induced by the formation of the slab. For the ZZ configuration obtained theoretically, we observed that the non-totally symmetric modes did not present any significant enhancement associated to the zigzag edge since their intensities are up to 3 orders of magnitude smaller than those of the totally symmetric modes. In contrast, mode A 1 g shows a clear Raman intensity enhancement due to the presence of the edge. A close view of mode A 2 g also shows a slight enhancement at the edge. For the XX configuration, there are no significant zigzag edge effects associated to the A 2 g mode. On the other hand, we clearly observed an edge enhancement for modes B 1 3g and B 2g . Finally, the XZ configuration shows that modes B 1g and B 2g possess the dominant intensities. In particular, one can clearly see a robust enhancement of the B 1g intensity at the edge. In addition, the totally-symmetric modes present weaker, albeit noticeable, enhancement effects at the zigzag edge.
In Supplementary Figure 8 we show the Raman intensities for the armchair edge with the ZZ, XX, and ZX configurations. For the ZZ configuration a noticeable intensity enhancement is observed for the totally-symmetric mode A 2 g and for the non-totally symmetric modes B 1g and B 2g . In addition, for configuration ZX the most enhanced mode is B 1 3g whereas mode B 2g does not show significant enhancement due to the edge.