Large area chemical vapour deposition grown transition metal dichalcogenide monolayers automatically characterized through photoluminescence imaging

Chemical vapour deposition (CVD) growth is capable of producing multiple single-crystal islands of atomically thin transition metal dichalcogenides (TMDs) over large areas. Subsequent merging of perfectly epitaxial domains can lead to single-crystal monolayer sheets, a step towards scalable production of high quality TMDs. For CVD growth to be effectively harnessed for such production it is necessary to be able to rapidly assess the quality of material across entire large area substrates. To date, characterisation has been limited to sub-0.1-mm2 areas, where the properties measured are not necessarily representative of an entire sample. Here, we apply photoluminescence (PL) imaging and computer vision techniques to create an automated analysis for large area samples of monolayer TMDs, measuring the properties of island size, density of islands, relative PL intensity and homogeneity, and orientation of triangular domains. The analysis is applied to ×20 magnification optical microscopy images that completely map samples of WSe2 on hBN, 5.0 mm × 5.0 mm in size, and MoSe2–WS2 on SiO2/Si, 11.2 mm × 5.8 mm in size. Two prevailing orientations of epitaxial growth were observed in WSe2 grown on hBN and four predominant orientations were observed in MoSe2, initially grown on c-plane sapphire. The proposed analysis will greatly reduce the time needed to study freshly synthesised material over large area substrates and provide feedback to optimise growth conditions, advancing techniques to produce high quality TMD monolayer sheets for commercial applications.


INTRODUCTION
Van der Waals layered crystals such as TMDs, hexagonal boron nitride (hBN) and graphene are defined by their strong covalent bonds in-plane and weak interlayer forces 1 . These characteristics allow for individual atomically thin layers to be easily removed from the bulk crystal, and for single layers to be brought together to build vertical heterostructures that display promising new properties [1][2][3] . As TMDs are exfoliated down to monolayer thickness their electronic band structure is altered due to dimensional confinement, causing the materials to become direct bandgap semiconductors that display photoluminescence (PL) under optical excitation 4,5 . MoSe 2 , WSe 2 , MoS 2 and WS 2 are documented as the most promising and widely studied TMD monolayers, emitting bright PL while remaining stable in air, at room temperature 6 . Such traits lend these materials to applications in optoelectronic devices 7 such as LEDs 3,8 , photovoltaic cells 9,10 , photodetectors 11 and single photo-emitters 12,13 . Their optical properties have attracted further attention due to the possibility for straightforward coupling to spin and valley degrees of freedom 14 . The extreme thinness and mechanical stability of these materials provide the potential for flexible and transparent devices 15 . However, to realise such devices, the goal of scalable and controlled production of mono-and multilayers must first be achieved, an objective which has been of great interest for commercial application ever since TMDs were rediscovered as atomically thin materials 4,5 .
The original method of monolayer production, micromechanical exfoliation 2 , produces high quality two-dimensional (2D) flakes, but with a low yield and random nature it has an inherent inability to be scaled up. This has lead to a catalogue of methods being developed including other top-down approaches such as liquidphase exfoliation 16 and various forms of thinning [17][18][19] , as well as bottom-up routes such as wet-chemical synthesis 20 and physical vapour deposition 21,22 . A further technique, chemical vapour deposition (CVD), has the capability to grow multiple monolayer islands across large area substrates and has been earmarked as the process that will deliver scalable production 23 .
Continuous lateral growth of these CVD islands leads to the inplane merging of single-crystal monolayers and eventually the formation of monolayer sheets. A problem encountered is that if randomly orientated, or irregularly shaped islands are allowed to coalesce, grain boundaries are formed at every merger of unaligned edges, which may negatively affect the optical and electrical characteristics of the sheet 24,25 . The density of grain boundaries can be reduced through epitaxial growth of the TMD islands on substrates with hexagonal symmetry, such as hBN [26][27][28] , mica 29 , c-plane sapphire 30,31 and graphene 32 . On these substrates, van der Waals interlayer interaction promotes the relative alignment despite the, often large, lattice mismatch 33 . In CVD the morphology of islands can be controlled to a degree by altering the local ratio of the transition metal and chalcogen atoms around a nucleation point. Typically the synthesised islands take on a hexagonal or equilateral triangle shape 34,35 . The regular and reproducible triangular morphology is of particular interest with respect to the production of continuous single-crystal films. On the majority of the crystalline substrates mentioned, hBN [26][27][28] , mica 29 and graphene 32 , two well defined groups of triangular islands rotated at 60°relative to one another are observed. The subsequent merging of these two groups results in either perfect stitching between objects at the same orientation, or the formation of mirror grain boundaries between those of opposite. On c-plane sapphire two more orientations were observed with a much lower probability, both at a 30°rotation from the first set 30 .
However, epitaxial growth is not a well-studied phenomenon for any TMD and substrate combination, and it is generally only assessed manually across sub-0.1-mm 2 areas, with analysis containing just hundreds of triangular crystals at most. Such small sample sizes are not necessarily representative of the CVD growth across an entire substrate, as all island properties, both physical and optical, vary with spatial position. If this ordered growth can be analysed across a larger area and a causality established with growth parameters, the gained feedback could be used for production of monolayer sheets with minimal grain boundaries, potentially realising commercial production of devices. Further requirement for regular orientation of islands arises from potential uses in van der Waals heterostructures, where a relative twist between adjacent layers is being actively researched as a new degree of freedom [36][37][38][39][40][41][42] .
Although the mechanisms that bring about the physical properties of CVD grown material are believed to be understood 20 , to date, there are no documented methods to provide feedback on an entire sample for optimisation of the desired properties such as, island size and density, quality of PL emission and orientation of individual islands. If the quality of CVD grown TMD monolayers is to improve and high-grade sheets and van der Waals heterostructures produced, it is first necessary to be able to reliably characterise entire large area samples.
Computer vision has been used in the field of 2D materials, predominantly as a tool for automatic mono-and multilayer identification of exfoliated material across large areas 43,44freeing up time taken searching for flakes. The technique of PL imaging has been developed in parallel, as a fast and innovative type of optical imaging, that can be simply applied to a standard optical microscope for wide-field and fast capture of fluorescent monolayer material 45 .
Here, these two techniques have been applied to create an automated analysis capable of characterising substrates on which CVD synthesis has been carried out. The size, orientation, relative intensity and homogeneity of PL emission are analysed using image processing software for each TMD island on substrates up to 65 mm 2 in size, containing more than 100,000 individual islands of one or two TMD materials. The analysis is demonstrated in application to two samples of WSe 2 grown directly onto hBN (sample 1 and 2, both 5.0 mm × 5.0 mm) and a third sample (sample 3) containing MoSe 2 and WS 2 islands on SiO 2 /Si (11.2 mm × 5.8 mm), where the WS 2 was grown directly onto the substrate while the MoSe 2 was grown on c-plane sapphire and subsequently transferred. The reported characterisation and analysis can be extended to any combination of TMDs and substrate where PL is emitted at room temperature. For WSe 2 monolayers grown directly on hBN, we observe two equally preferential orientations, while MoSe 2 samples originally grown on c-plane sapphire and subsequently transferred onto SiO 2 /Si substrate demonstrate four prevalent preferential orientations.

RESULTS
Imaging and image processing framework We use computer vision functions from the MATLAB image processing toolbox 46 to analyse the high-resolution PL image of the entire sample on a large substrate of a few 10 s of mm 2 , acquired by manually stitching the images taken with a ×20 magnification microscope objective. The monolayer TMD islands were identified through their PL as the only luminescent objects on the studied substrates, and their size and shape were extracted using MATLAB shape recognition functions. Only monolayer islands with the shape close to equilateral triangle were used for orientation analysis, while monolayer islands with all morphologies were used for size, mean PL intensity, and material coverage analysis, as described below in this section, and further detailed in the 'Methods' section and in Supplementary information. Figure 1a shows a schematic diagram of how a dichroic mirror can be used to isolate PL emission in a generic microscope set-up, reflecting the visible light while transmitting in the near infrared 45 Fig. 1 Imaging and image processing of CVD grown TMD monolayers. a A schematic showing how PL emission is isolated from the collected light in a typical BF microscope set-up using a 550-nm short-pass filter, 550-nm long-pass dichoric mirror and 600-nm long-pass filter combination. b-d A comparison between images at ×50 magnification obtained using BF (b), DF (c) and PL (d). In (d), the two materials are outlined in light blue with MoSe 2 appearing as a light pink colour and WS 2 as a dark red. The image was taken with a 10 s exposure time and 9.6x analogue gain. The red scale bars in each image are 50 μm. The binary image resulting from the image processing of (d) is shown in (e), with each individual object highlighted by their central point (centroid) shown by a blue marker, and their extrema shown by red dots. f Illustration of how the applied image analysis programme attempts to find a triangular outline of a single isolated island by combining or discarding the extrema, resulting in three clear combined points, E1-3. See further details in text and Supplementary Information (Supplementary Note 1). g A visual representation of how the programme calculates the orientation of each island from the positions of the combined extrema, E1-3. The light blue equilateral triangle is fitted to each combinations of two points and the angle, θ, along with an uncertainty, is calculated.
(see details in the 'Methods' section). A comparison between images of monolayers at ×50 magnification in bright-field (BF), dark-field (DF) and PL modes is shown in Fig. 1b, c and d, respectively, for the same area of the MoSe 2 -WS 2 sample on SiO 2 / Si (referred to as sample 3 below). In the BF image objects are distinguishable from the purple substrate, although contrast is low and the two materials are difficult to discriminate. In the DF image edges are well defined, but other features on the sample surface, such as dust, appear bright-masking monolayers. On top of this, the edges of both materials appear as the same colour and objects of WS 2 are only recognisable from those of MoSe 2 by their size under careful analysis. The PL image shows the greatest contrast between the substrate and materials as only the latter have the potential to display PL, the presence of which confirms their monolayer thickness. The material dependent wavelength of the emitted light at room temperature allows multiple TMDs to be easily identified in a single image and analysed separately 45 , with PL peaked at 785 nm for MoSe 2 , 750 nm for WSe 2 , 675 nm for MoS 2 and 630 nm for WS 2 , where wavelengths above 700 nm appear as false colours. The boundaries of all isolated islands are clear and the dominant triangular morphology is evident 35 . Compared with DF and, to a degree, BF, unwanted contaminants such as dust particles have little effect on the quality of a PL image.
In order to measure properties of the islands, size, orientation, density and relative PL emission, it is necessary to first convert the PL image into a binary form, with pixel values set to either on (1), representing PL emission and thus a monolayer, or off (0), for areas of substrate or thicker material. This is achieved through colour thresholding to isolate all red and light pink coloured objects in an image, as shown in Fig. 1e. In this form, image processing functions, taken from MATLAB's image processing toolbox 46 , can be applied to improve the quality of analysis (see details in Supplementary Note I). The chosen parameters of the feature detection algorithm employed for the images reported in this work allow measurement of islands with the lateral size roughly above 6 μm. The full colour 8-bit image is also converted to grayscale, where each pixel is represented by a single value between 0 and 255. This pixel value is calculated from the channels of the given format that define its colour before the conversion, in this case red, green and blue (RGB). The properties of islands are measured from both binary and grayscale images. Size and extrema, which are used to find orientation, are taken solely from the binary image, and mean pixel value and standard deviation (STD) of pixel value are extracted from the masked grayscale image. The latter two properties provide a relative measure of the average brightness and homogeneity of PL emission across an objects surface.
To measure the orientation of each object the assumption is made that all islands are equilateral triangles, the most common shape seen in Fig. 1b-e. As mentioned in the introduction, this is a well documented morphology whose growth mechanism and parameters are reasonably well understood 35,47,48 . The positions of the extrema, shown in red in Fig. 1e, f, are called from the 'regionprops' MATLAB function, returning the furthest eight points on an objects surface from its centroid. These are averaged together based on proximity to one another and size of the given island. Only the vertexes that are defined by at least two of the extrema returned by the function are considered for the analysis, thus ensuring an unambiguous identification. If an island is found to have anything but three vertexes it is not a triangle and is discarded. Figure 1f demonstrates this step, where the positions of the extrema are marked by red stars and the position of the centroid is marked by a blue asterisk. Those extrema contained in a blue oval would be combined to an average position, and the one circled in red would be discarded. An equilateral triangle is then fitted to each combination of two vertexes and an orientation found. Figure 1g shows this step, where the light blue equilateral triangle is fitted to the points E1 and E2, and θ is the calculated orientation. Typically an island grown via CVD is not perfectly equilateral, and so fitting an equilateral triangle to each side gives us a range of orientations which can be used to find an uncertainty. This, in turn, can be applied to filter out islands displaying incorrect morphology by discarding any measurement with a percentage uncertainty >30% of the mean value for the data set. A more detailed explanation of the image processing and the steps of the analysis can be found in the Supplementary Note I and in Supplementary Fig. 1. Figure 1d and e shows some limitations of the analysis. If two single-domain islands grow into one another the proposed analysis fails to identify the position of the grain boundaries from the PL image. Thus, in the binary image, the two islands will be treated as a single object whose orientation data will be discarded. However, across an entire substrate of dimensions 11.2 mm × 5.8 mm, seen in Fig. 2a, an orientation could be measured for 28.2% of the objects-proving to be sufficient to gauge the level and strictness of aligned growth. It is important to note that only the orientation analysis discards data, all other properties such as size, mean pixel value, STD of pixel value and island density, are taken from all monolayer area, regardless of shape. This orientation measurement method could be easily adapted for a hexagonal morphology.
Large area characterization of a MoSe 2 -WS 2 sample A macroscopic image of the MoSe 2 -WS 2 sample on SiO 2 /Si can be seen in Fig. 2a, showing the presence of material but no details of individual islands. The entire sample was manually mapped in PL imaging mode at ×20 magnification with 10 s exposure time and 9.6x analogue gain, as described in Methods, and the set of image processing tools outlined above were applied to each frame (373 μm × 497 μm). Details of the fabrication of this sample can also be seen in Methods.
Heatmaps displaying the number of isolated MoSe 2 and WS 2 objects contained in each image, represented by a single tile on the map, are shown in Fig. 2b and c. The size of all islands in each image are averaged and represented by a tile on the heatmaps of Fig. 2d and e. It can be seen that MoSe 2 islands are much more homogeneous in size across the surface of the substrate, although flakes do get larger closer to the middle of the left edge. WS 2 vary more in size, and the areas containing the largest islands appear to roughly align with the sections of highest number density. What is quite clear from the heatmaps, and in agreement with the reference image Fig. 2f, whose position on the maps is outlined by a black box, is that MoSe 2 islands are small and densely packed together while the WS 2 islands are much larger and have a lower nucleation density. MoSe 2 is seen across two-thirds of the substrate from left to right while WS 2 has mostly developed along the top.
The histogram, Fig. 2g, displays a distribution of size with a logarithmic scale for all objects on the sample. This shows that the WS 2 tails off at much larger areas and that the smallest islands are dominant in number for both materials. Both distributions are offset from 0 since objects, <200 pixels (19.3 μm 2 ) in size for MoSe 2 and 300 pixels (28.9 μm 2 ) for WS 2 , were omitted in the image analysis.
If a monolayer sheet completely covering the substrate is considered 100% coverage, the identified 97,672 isolated MoSe 2 objects cover 6.90% of the substrate and 2573 WS 2 objects cover 0.77%. If the 39.4% of the substrate that is completely empty of any material, shown by the white area at the bottom right of the heatmaps, is omitted from the calculation, then the total coverage is 11.0% and 1.24% for MoSe 2 and WS 2 , respectively. These coverage values are slightly underestimated as they do not include small objects removed during the image processing.
The optical properties of a sample can be displayed in much the same way as the size and density information, to reveal how PL emission changes across the surface of a sample. The heatmaps in Fig. 3 show how the average of the mean pixel value for each object in an image changes with position on the sample for MoSe 2 (a) and WS 2 (b). In the PL imaging mode the mean pixel value for an object is directly related to the strength of that object PL emission, the camera settings used to capture the image and the optical excitation power.
Relative comparison between the properties of individual material within a heatmap is possible, as all images were acquired with identical gain and exposure settings, and with the same excitation power. For MoSe 2 , Fig. 3a shows clear gradients between areas of the highest and lowest PL emission, while WS 2 in (b) shows a more random nature with respect to spatial position. The distributions seen in Fig. 3c can shed some light on the variation in the average PL emission of individual islands across the sample through their full width half maximum (FWHM), measured as 28 for both materials. The standard deviation (STD) from the mean pixel value of an object tells us how much these values vary across its surface and therefore how homogeneous its PL emission is. As a distribution, Fig. 3d, STD allows the consistency of the optical emission to be seen on an island by island basis, across an entire substrate. Here, the PL emission homogeneity of WS 2 islands may be slightly compromised by the quenching that occurs where heterostructures are formed with the smaller islands of MoSe 2 45 . The PL quenching is caused by the interplay of efficient charge and energy transfers between the two materials 49,50 and fast non-radiative decay processes 51,52 . For this reason, in the current form, the algorithm does not recognize the MoSe 2 islands overlapping with the WS 2 monolayers. The number of such MoSe 2 islands is statistically insignificant and does not impact the conclusions of the characterization analysis for our sample.
Measurement of epitaxial growth degree The strictness of epitaxial growth to the hexagonal lattices of their substrates has been assessed for the three samples introduced earlier. Samples 1 and 2 are of WSe 2 grown directly onto hBN, both 5.0 mm × 5.0 mm in size. Sample 3 is of MoSe 2 -WS 2 on SiO 2 / Si, 11.2 mm × 5.8 mm in size, where MoSe 2 was grown on c-plane sapphire and then transferred, and the WS 2 was grown directly onto the SiO 2 /Si. The details of how these samples were grown can be found in Methods. The samples of WSe 2 were manually mapped in PL imaging mode at ×20 magnification, in a similar way as described above for sample 3. The image processing tools were applied to each frame and the shape of every identified object was analysed, and an orientation measured, when possible, via the method introduced earlier. The histograms in Fig. 4b, d and f display the results of these measurements.
We further visualise and support our findings on the singlecrystal domain orientation presented in Fig. 4b, d and f by using Fast Fourier transforms (FFTs) analysis of the microscopy images. An FFT represents an image as a collection of two-dimensional sinusoidal waves of varying wavelength. A commonly orientated boundary, such as one side of a repeated shape in an epitaxial sample, would be represented by a line running at an angle perpendicular to the boundary edge through the centre of the FFT, assuming no pattern in nucleation position or size. For a triangular morphology and two dominant directions of alignment we thus would expect three lines passing through the centre of the image, whereas for four dominant orientations, six lines should be observed. Each of the FFT graphs shown in Fig. 4 represent the entirety of each sample. These images were created An example image of sample 1 can be seen in Fig. 4a, and on inspection no preferred orientation of individual islands can be seen. This is in agreement with the FFT, inset to Fig. 4a, where there are no obvious lines observed, instead a circular shape is seen which is characteristic of randomly orientated islands. Across sample 1, an orientation could be measured for 1419 objects from a possible 5272 that were found by the programme. Thus, 26.9% of objects could be measured which is 5.65% of all monolayer material by area. The low percentage of monolayer area where orientation is measurable is due to many different types of growth observed in this sample, including multi-pointed stars 35 , more irregular shapes and monolayer sheets, all composed of multiple crystal domains, which tend to be much larger objects than a single triangular island. The image shown in Fig. 4a is an example of the areas which can be most effectively analysed. The histogram of orientation in Fig. 4b shows a bi-modal distribution with peaks centred around 35°and 90°, which is close to the 60°s plit expected from a sample with two directions of alignment at 180°to one another and a triangular morphology with three-fold symmetry. The two distributions are very broad suggesting only a small degree of alignment to the hBN substrate. The orientation distribution expected from typical epitaxial growth measured for sample 2 of WSe 2 on hBN is seen clearly in Fig. 4d. An example image of this highly aligned sample is seen in (c), along with an inset FFT. The orientation distribution in Fig. 4d shows two peaks centred at 37°and 95°, with an almost equal number of islands contained in each. The degeneracy in number of islands at each orientation is a well documented characteristic of epitaxial growth 27,28,30 . The narrowness of the peaks display how closely aligned the islands are to the crystal lattice of hBN with FWHMs of 9°and 7°for the 37°and 95°peaks, respectively. The FFT shows three clear lines as expected. This sample was identical in size to sample 1, and 5576 isolated single-domain islands were measured from a possible 14,033 objects, 39.7% by number but only 9.90% by area. Similarly to sample 1, sample 2 contains large monolayer sheets, present on nearly a third of the sample, and irregularly shaped islands whose orientation cannot be measured through this method. However, triangular objects were a more consistent equilateral shape than in sample 1, making them easier to identify and measure. The image seen in Fig. 4c shows an area where the orientation can be seen clearly, but it is not representative of the entire sample.
The final orientation histogram, Fig. 4f shows the same analysis applied to sample 3 of MoSe 2 -WS 2 described in Figs 2 and 3. The WS 2 islands are low in number and only 263 of the 2573 identified could be measured, 10.2% of objects by number and 11.1% by area. This small number of objects being measured is due to many of the isolated islands not being triangular, and of those that are, many have their PL quenched by the formation of heterostructures with MoSe 2 -affecting their observed shape. The flat distribution suggests a lack of epitaxial growth for this material, although the relatively small sample size means any conclusions drawn are unlikely to be reliable. The MoSe 2 islands show a more ordered behaviour. Of the 97,672 isolated MoSe 2 objects identified, 27,516 could be measured by the orientation analysis programme-28.2% of objects by number and 23.2% of the material by area. The much higher percentage of monolayer area that could be measured, compared with previous two samples, is due to this sample containing no areas of continuous monolayer sheets. However, irregular shapes from merged domains are still common. The image in Fig. 4e is a good example of the average growth across the covered sections of the sample. By eye the MoSe 2 islands appear to be randomly orientated, but the distribution shows four peaks spaced roughly 30°apart at 12.5°, 42.0°, 73.0°and 103°. These results can be reconciled with the synthesis of the materials, as WS 2 was grown on SiO 2 /Si and MoSe 2 on the hexagonal crystal structure of c-plane sapphire. These four preferential directions have been observed previously on c-plane sapphire across a 220 μm × 220 μm area 30 , where two of the orientations-60°apart-were dominant, shown by 91.5% of the islands, while the other two, 30°from the first were only seen in 6%. Here, although we see a more dominant peak at 12.5°, the four clear peaks highlight the existence of four significant preferential directions of growth. It appears that the strictness of alignment is lower for this sample than for sample 2 and FWHMs cannot be found accurately, but the pattern can still be observed in the FFT, inset to Fig. 4e, showing six lines, representing the commonly orientated boundaries of the material. For the CVD growth of monolayer MoS 2 on mica 29 , WS 2 and MoS 2 on hBN [26][27][28] , and WS 2 on graphene, two preferential directions have been seen in previous literature with triangular islands rotated at 60°relative to one another. Both two 31 and four 30 orientations have been observed for growth on sapphire. For the latter, the second set, 30°from the first set, had been shown to have a low probability 30 . To further reduce grain boundaries, it has been discovered that the introduction of atom vacancies in hBN acts to trap W atoms of WSe 2 , breaking the typical degeneracy in the number of islands at either orientation 53 . However, assessment of aligned growth has so far been limited to tens or hundreds of islands across sub-0.1-mm 2 areas, and therefore conclusions of large area epitaxial arrangement are hard to make considering the relatively unpredictable nature of CVD growth. Although previous work demonstrates that some degree of control is possible in CVD growth, having a means to assess the strictness of the alignment across a whole sample would allow growers to harness this control to reduce the density of grain boundaries in monolayer sheets, potentially improving electrical and optical properties of TMD films. Using image processing, we measure the strictness of epitaxial growth across a 5.0 mm × 5.0 mm WSe 2 /hBN sample, where alignment was obvious on inspection, and uncover four prevalent orientations in MoSe 2 initially grown on c-plane sapphire, across a 11.2 mm × 5.8 mm MoSe 2 -WS 2 sample on SiO 2 /Si.
There are limitations which must be noted. The merging of single-domain islands causes data on epitaxially grown material to be discarded, along with a possible overestimation of size measurements and an underestimation of nucleation density. The effects on size and density can be countered with a further image processing technique called segmentation, although the artificial borders created during this process are not always accurate and it tends to overcompensate for the problem. On the other hand, for orientation, in most cases a large enough number of objects can be detected across a substrate to give a representative measurement of the level of epitaxial growth. Any technique that relies on PL imaging will encounter similar problems and a more complex method would be necessary to extract orientation and grain boundaries from merged domains. As can be seen in the Supplementary Fig. 2, the degree of epitaxial growth of objects can still be assessed in samples of 2D materials that do not emit PL at room temperature using DF imaging, providing the correct morphology is present.

DISCUSSION
A new automated analysis has been developed capable of measuring the characteristics of monolayer islands across tens of mm 2 areas of CVD grown TMD semiconductors, from images that map the sample surface. The analysis uses PL imaging techniques and can be applied to various combinations of semiconducting TMD monolayer and substrate which emit PL at room temperature. Functions from MATLAB's image processing toolbox are used to measure island size and density, as well as relative average value, and homogeneity of, PL emission. Importantly, for each isolated monolayer island across a sample we measure orientation, and obtain data on the presence and degree of epitaxial growth. An FFT analysis has been further applied to uncover any underlying patterns in the sample and to verify the results of the orientation histogram. Four preferential directions of alignment, not apparent to the naked eye, were found in MoSe 2 islands grown on c-plane sapphire, a characteristic that had not been seen previously to such a degree. This analysis method coupled with the automated mapping of CVD samples in PL image mode forms a fast and reliable characterisation tool.
With the use of dark-field imaging the analysis of epitaxial growth can be extended to many other materials that do not emit PL. Such analysis will reduce the amount of time needed to study CVD samples or 2D materials produced by other scalable techniques, and has particular application to developing a consistent method of production for high quality 2D monolayer sheets and heterostructures. Applications in other thin film material systems, such as perovskite or organic semiconductors are also possible.
Image analysis will likely become the norm as progression in the field leads to inevitable commercial production of devices that require a consistent quality of TMD monolayer, without time for constant manual inspection. The utility of the automated analysis developed in this work comes in its ability to characterise the quality of grown material automatically and unearth patterns that would be missed by manual inspection, such as the four preferential directions found in material grown on c-plane sapphire (Fig. 4f), and the quantitative level of homogeneity in MoSe 2 island size across the same sample (Fig. 2d).
With thousands of objects analysed, a statistical approach to assessing the quality of an entire substrate can be taken, as well as the ability to focus on the properties of individual islands such as, size, orientation and quality of PL emission. The analysis holds vast amounts of information, providing a deep understanding of the characteristics of a sample and the opportunity to quantitatively compare samples grown under different conditions.

METHODS
Growth of single-layer WS 2 WS 2 was grown on a SiO 2 substrate via CVD. Before the growth of WS 2 , SiO 2 was cleaned by a piranha solution and was spin-coated by sodium cholate which acted as a seeding promoter. WO 3 (99.998%, Alfa Aesar) was dissolved in the water/ammonia solution (9:1) and 1 mL of the solution was coated on the crucible. This crucible was placed at the centre of the furnace and 100 mg of S (99.999%, Alfa Aesar) was placed at the upstream entry of the furnace. Then, the temperature of the tube furnace was increased up to 900°C for 24 min under a steady flow of Ar gas (100 sccm) in the ambient condition. When the furnace reached 600°C the S vaporised. Then, temperature of the tube furnace was maintained at 900°C for 30 min for the WS 2 growth. Afterwards, the tube furnace was cooled down to room temperature under the Ar flow.
Growth of single-layer MoSe 2 and transfer to SiO 2 /Si MoSe 2 was grown on c-plane sapphire by CVD. Two precursors, MoO 3 (99.97%, Sigma-Aldrich) and Se (99.999%, Alfa Aesar), were used for the growth. One hundred and fifty milligrams of Se was placed at the upstream entry of the furnace and 60 mg of MoO 3 powder was placed at the centre of the furnace. A crucible containing MoO 3 was partially covered by a SiO 2 /Si wafer to reduce intense evaporation of the precursor. The sapphire substrate was located next to the crucible that contained MoO 3 . Before the tube furnace was heated, the tube was evacuated for 30 min and filled with the Ar gas achieving ambient pressure. The temperature of the furnace was increased up to 600°C for 18 min under a steady flow of Ar gas (60 sccm) and H 2 gas (12 sccm). When the furnace reached 600°C, Se was vaporised by heating the upstream entry of the tube up to 270°C using a heating belt. Finally, temperature of the tube furnace was increased to 700°C and maintained for 1 h for the MoSe 2 growth. Afterwards, the tube furnace was cooled down to room temperature while the Ar flow was maintained without H 2 . To transfer MoSe 2 on top of the SiO 2 /Si substrate containing WS 2 , polystyrene was used to maintain the sample quality instead of poly(methyl methacrylate) that has been widely used previously 54 .

Growth of single-layer WSe 2 on hBN
To fabricate WSe 2 on hBN, multilayer hBN was initially grown on c-plane sapphire 55 . WO 3 (99.998%, Alfa Aesar) and Se (99.999%, Alfa Aesar), were used for the WSe 2 growth. Three hundred milligrams of Se was placed at the upstream entry of the furnace and 120 mg of WO 3 powder was placed at the centre of the furnace. To reduce the influence of humidity, a small amount of NaCl was mixed with WO 3 powder. The multilayer hBN on sapphire substrate was positioned next to the crucible containing WO 3 . Before the tube furnace was heated, the tube was evacuated for more than 30 min. Then, the temperature of the tube furnace was increased to 800°C for 24 min under a steady flow of Ar gas (120 sccm) and H 2 gas (20 sccm). When the furnace reached 800°C, the Se was vaporised by heating the upstream entry of the tube to 270°C using a heating belt. Finally, temperature of the tube furnace was increased to 870°C and maintained for 1 h for the WSe 2 growth. Afterwards, the tube furnace was cooled down to room temperature under Ar flow.

PL imaging
The PL images analysed, such as Fig. 1d, were taken with an adapted industrial microscope (LV150N, Nikon). A 550-nm short-pass filter is positioned in the path of the white light from the illumination source and a 600-nm long-pass filter is placed before the camera (DS-Vi1, Nikon) to isolate PL emission. A further 550 long-pass dichroic mirror is applied to direct the excitation light and collected PL emission. A detailed description of the method can be found in ref. 45 . Using this technique, the samples discussed in reference to Fig. 1 through 4 were mapped manually at ×20 magnification. The images of WSe 2 on hBN were taken with 4 s exposure time and an analogue gain of 5.6x, those of MoSe 2 -WS 2 on SiO 2 /Si were taken with 10 s and 9.6x, respectively. Each image was 373 μm × 497 μm in area.

Image processing
Images were analysed in MATLAB using functions from the image processing toolbox 46 . The colour thresholding application was used to isolate monolayer material in a PL image and was applied to each combination of monolayer and substrate. The 'regionprops' function was used for the analysis of size, island density and pixel value with minimal custom code. The programme that measured the orientation of equilateral triangles was developed specifically for the analysis of epitaxial growth. It took the programme 43 s to analyse the orientation of the images for sample 1 containing 5272 objects, 69 s to analyse that of sample 2 of 14,033 objects, and 1 h 8 min and 13 s to analyse sample 3 containing 100,245 objects of two different colours. FFTs, representing entire substrates, were created by superimposing individual image FFTs and were subsequently artificially enhanced for clarity by transforming the pixel values using contrast-limited adaptive histogram equalisation. Further details of the image processing and a more complete explanation of the analysis can be found in the Supplementary Note I.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.