Lorenzo Valzania
email: lorenzo.valzania@lkb.ens.fr

Jonathan Dong

Sylvain Gigan

Accelerating ptychographic reconstructions using spectral initializations

Ptychography is a promising phase retrieval technique for label-free quantitative phase imaging. Recent advances in phase retrieval algorithms witnessed the development of spectral methods, in order to accelerate gradient descent algorithms. Using spectral initializations on experimental data, for the first time we report three times faster ptychographic reconstructions than with a standard gradient descent algorithm and improved resilience to noise. Coming at no additional computational cost compared to gradient-descent-based algorithms, spectral methods have the potential to be implemented in large-scale iterative ptychographic algorithms.

I. INTRODUCTION

Ptychography is a computational imaging technique that enables label-free, quantitative phase imaging [START_REF] Rodenburg | Ptychography and related diffractive imaging methods[END_REF]. It is based on a simple principle: scan a probe across a sample, collect the corresponding intensity diffraction patterns (also known as 'ptychograms'), and reconstruct an image of the object of interest. Because it does not require complex optical elements, it has been adapted to a variety of settings and spectral ranges, from electron microscopy [START_REF] Jiang | Electron ptychography of 2D materials to deep subångström resolution[END_REF], for which it was originally conceived in the late 60s [START_REF] Hoppe | Diffraction in inhomogeneous primary wave fields. 1. Principle of phase determination from electron diffraction interference[END_REF], to X-rays [START_REF] Pfeiffer | X-ray ptychography[END_REF] and extreme ultra-violet light [START_REF] Zhang | Tabletop nanometer extreme ultraviolet imaging in an extended reflection mode using coherent Fresnel ptychography[END_REF], all the way down to the terahertz (THz) range [START_REF] Valzania | Terahertz ptychography[END_REF].

The computational reconstruction in ptychography requires to solve a phase retrieval problem, where the phase of the diffracted electric field has to be recovered from intensity-only measurements. Such a problem is tractable if a minimum overlap ratio (empirically estimated around 60% [START_REF] Bunk | Influence of the overlap parameter on the convergence of the ptychographical iterative engine[END_REF]) is ensured between subsequent probe positions. However, as a non-linear and nonconvex optimization problem, it is still not completely understood and convergence towards its global minimum is not guaranteed.

One way to avoid local minima is to provide an initial estimate already close to the solution, and further refine it with iterative algorithms. To this end, spectral methods have recently been proposed for the general phase retrieval problem with independent and identically distributed (i.i.d.) random measurements, solved with a gradient descent (GD) approach [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF]. The initial estimate is defined as the leading eigenvector of a covariance matrix, constructed from the experimental measurements and the acquisition parameters. The quest for improved reconstruction performance culminated in the derivation of optimal spectral methods for the random measurements setting [START_REF] Mondelli | Fundamental limits of weak recovery with applications to phase retrieval[END_REF][START_REF] Luo | Optimal spectral initialization for signal recovery with applications to phase retrieval[END_REF]. Due to these breakthroughs, spectral methods for phase retrieval are rather well understood theoretically.

Although spectral initializations to solve ptychography have been showcased on simulated measurements [START_REF] Marchesini | Alternating projection, ptychographic imaging and phase synchronization[END_REF], no gain was ever reported when employed on experimental data. Here, we provide the first experimental demonstration of ptychographic reconstructions improved and accelerated by spectral initializations.

II. METHODS

Let us begin by describing a ptychogram y (l) (x) with the following forward model

y (l) (x) = |P d {a(x -x (l) ) ψ(x)}| 2 , (1)
where x is a two-dimensional spatial coordinate, ψ(x) is the complex transmission function of the object, scanned with the probe a(x) at the positions x (l) for l = 1, . . . , L, with L being the total number of images. P d is a known linear operator describing the transmission through the optical system of optical length d between the object plane and the detector plane. We hereby point out that the model equally applies to Fourier ptychography [START_REF] Konda | Fourier ptychography: current applications and future promises[END_REF], if one regards ψ(x) and a(x) as the Fourier transform of the object transmission and the probe functions respectively, and P d is an inverse Fourier transform. In what follows the case of ptychography will be considered, with a setup sketched in the top row of Fig. 1.

To recover the object, ptychography solves a phase retrieval problem that can be formulated in the vectorized form

y = |Sψ| 2 .
(2)

Leading eigenvector of 𝒁 Iterative optimization (GD, PIE, etc.)

FIG. 1. Sketch of a ptychography setup and algorithmic pipeline with our spectral method. Top: an aperture generates the probe a(x), that is scanned with shifts x (l) across an object of interest ψ(x), for l = 1, . . . , L. For each probe position, a camera records a ptychogram y (l) (x) after free space propagation by a distance d. Middle: the ptychograms are preprocessed and used to construct a weighted covariance matrix Z. Bottom: the spectral estimate, obtained as the leading eigenvector of Z, is already close to the solution and is further refined using an iterative optimization (dark blue point). Initial estimates further away from the solution (light blue point) may be stuck in local minima. GD: gradient descent; PIE: ptychographic iterative engine.

Here, ψ ∈ C n represents the vector of the unknown object transmission coefficients at each of its n pixels, S ∈ C p×n is the sensing matrix, with p being the number of measured intensity values (i.e., p = Lm, where m is the number of camera pixels), y ∈ R p is the vector collecting the ptychograms, and | • | denotes the elementwise modulus operation. When the entries of S are i.i.d. random variables, ψ is initially estimated as the leading eigenvector of the following n × n matrix (Fig. 1, middle row):

Z ≡ S † diag{T (y)}S , (3)
where S † denotes the Hermitian conjugate of the matrix S, T is a preprocessing function acting element-wise on the measurements, and diag{u} denotes a diagonal matrix with a general vector u on its main diagonal. A widely used method to compute the leading eigenvector is power iterations, that converge exponentially towards the eigenvector corresponding to the largest eigenvalue in absolute value [START_REF] Dongarra | Numerical linear algebra for highperformance computers[END_REF]. This estimate is finally used as the first estimate of an iterative optimization algorithm like gradient descent or, in the case of ptychography, dedicated routines like the ptychographic iterative engine (PIE) [START_REF] Rodenburg | Ptychography and related diffractive imaging methods[END_REF]. The reconstruction pipeline is graphically summarized in the bottom row of Fig. 1.

In order to apply spectral methods to ptychography, the forward model in the continuous domain of Eq. (1) needs to be brought into the discrete vectorized form of Eq. ( 2). Let a (l) ∈ C m×n be the discrete and vectorized form of a(x -x (l) ), and P d ∈ C m×m be the matrix computing the linear transform P d . Now we can adopt the matrix factorization used in [START_REF] Marchesini | Alternating projection, ptychographic imaging and phase synchronization[END_REF] and finally write y = |Sψ| 2 , with S ≡ P A and

P ≡    P d • • • 0 . . . . . . . . . 0 • • • P d    ∈ C Lm×Lm , A ≡   
a (1) . . .

a (L)    ∈ C Lm×n , (4
) from which, using Eq. ( 3) and the unitarity of P (i.e., P † = P -1 ), one obtains

Z = A † P -1 diag{T (y)}P A . ( 5
)
Our choice of the preprocessing function, reading

T (y ) ≡ max(-β, 1 -1/y ) , (6)
where y collects the ptychograms, each normalized to its own mean intensity, and 0 < β < 1, was inspired by the theoretical works in [START_REF] Mondelli | Fundamental limits of weak recovery with applications to phase retrieval[END_REF][START_REF] Luo | Optimal spectral initialization for signal recovery with applications to phase retrieval[END_REF][START_REF] Wang | Solving systems of random quadratic equations via truncated amplitude flow[END_REF]. For y 1, the function corresponds to the optimal preprocessings of [START_REF] Mondelli | Fundamental limits of weak recovery with applications to phase retrieval[END_REF][START_REF] Luo | Optimal spectral initialization for signal recovery with applications to phase retrieval[END_REF]. The lower bound set at low intensity values prevents large negative eigenvalues of Z, that will degrade the performance of the power iteration algorithm. It is indeed known that power iterations return the eigenvector corresponding to the largest eigenvalue in absolute value. In addition, it reminds of the truncated amplitude flow algorithm, whose resilience to noise has been proved [START_REF] Wang | Solving systems of random quadratic equations via truncated amplitude flow[END_REF]. After plugging Eq. ( 6) into Eq. ( 3), the power method resembles a GD update, proving that no computational complexity is added when replacing one GD iteration with one power iteration (see Appendix B).

In the box Algorithm 1 we summarize the overall algorithm, consisting of the combination of a spectral initialization and a GD optimization with amplitude loss function

L(ψ) ≡ √ y -|Sψ| 2 2
, where • 2 indicates the L 2 -norm. The corresponding code is available at Ref. [START_REF] Dong | Ptychography with a spectral initialization[END_REF], and more details on a memory-efficient implementation of power iterations are outlined in Appendix C.

III. EXPERIMENTS

Experiments were performed using the THz imaging setup at Empa, the Swiss Federal Laboratories for Materials Science and Technology, equipped with a farinfrared gas laser (FIRL 100, Edinburgh Instruments, Livingston, Scotland) emitting several tens of mW of continuous-wave power at the wavelength λ = 96.5 µm. An uncooled microbolometer array detector featuring m = 480 × 640 pixels on a pitch of 17 µm (Gobi-640-GigE, Xenics, Leuven, Belgium), was used as a THz camera Algorithm 1: Solve ptychography with a spectral initialization

Input: Measurements y, sensing matrix S, preprocessing function T , initial estimate ψ0, number of power iterations M , number of GD iterations N , step size γ 1 Z = S † diag{T (y)}S /* Spectral initialization */ 2 for t = 1, . . . , M do

3 ψt = Zψt-1 4 ψt = ψt/ ψt /* Gradient descent */ 5 for t = M + 1, . . . , M + N do 6 ∇L(ψ) = S † Sψt-1 -S † diag{ √ y/|Sψt-1|}Sψt-1 7 ψt = ψt-1 -γ∇L(ψ)
8 return ψM+N [START_REF] Hack | Comparison of thermal detector arrays for off-axis THz holography and real-time THz imaging[END_REF]. The object was a 2-mm-thick polypropylene (PP, refractive index n P P = 1.51 and absorption coefficient α P P = 1.5 cm -1 at λ = 96.5 µm [START_REF] Lee | Principles of terahertz science and technology[END_REF]) slab where three intersecting rings were engraved by laser ablation at depths in the range 37-234 µm. Fig. 1 shows the wrapped phase shift distribution induced by the object at λ = 96.5 µm.

We let a plane wave diffract through a circular aperture with a diameter of 3 mm and propagate by 12 mm before impinging on the object. The object was scanned across a grid of 21 × 7 points at overlap ratios of 87% and 80% [START_REF] Bunk | Influence of the overlap parameter on the convergence of the ptychographical iterative engine[END_REF] along its vertical and horizontal axis in Fig. 1, respectively. The ptychograms were recorded after a free space propagation of d = 9.1 mm. This resulted in a Fresnel number around 2, making the angular spectrum propagator a suitable choice for computing P d [START_REF] Valzania | Terahertz ptychography[END_REF][START_REF] Goodman | Introduction to Fourier optics[END_REF].

IV. RESULTS

Although spectral methods in combination with preprocessing functions similar to Eq. ( 6) are well grounded in the theory of random sensing, their effectiveness and resilience to noise are also maintained when S is not random to a certain extent, as in our experimental setting. To prove this, we analyzed three datasets, collected at different signal-to-noise ratios (SNR) by duly tuning the laser power and denoised with the procedure borrowed from [START_REF] Williams | Fresnel coherent diffractive imaging: treatment and analysis of data[END_REF]. Results for each noise level are compiled in the columns of Fig. 2, labelled with their corresponding SNR, calculated as the ratio between the mean intensity in the center of each ptychogram and that at its outermost pixels, where the diffracted intensity is negligible. Starting from a flat estimate in both amplitude and phase, the power method yielded the spectral estimates in row (a) (the normalized reconstruction error [START_REF] Maiden | An improved ptychographical phase retrieval algorithm for diffractive imaging[END_REF][START_REF] Maiden | Further improvements to the ptychographical iterative engine[END_REF] is given above each image). Although we cannot guarantee that the leading eigenvector of Z is the sought solution even in the absence of randomness, we have indeed observed a decrease of the loss L in the early iterations of the power method, which can be explained by comparing one iteration of the power method with one GD iteration (see Appendix B). When such a decrease came to a halt, power iterations were stopped, and the obtained spectral estimate was used as the initial guess for GD, run with an exponentially decreasing step size, compatible with a backtracking line search, to ensure stable convergence. At all noise levels, two GD iterations are enough to provide a satisfactory reconstruction of the object, shown in row (b). For comparison, we also performed reconstructions with a standard Wirtinger flow GD algorithm without spectral initialization, using the same initial estimate provided to the power method in row (a) and the step size recommended in [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF]. The results, after the same total number of iterations run in row (b) (summing power iterations and GD iterations), are shown in row (c). Note that in the noisiest case, the reconstruction is dominated by high frequency artifacts at the periodicity of the acquisition scan, caused by a signal bias in the measurements [START_REF] Thurman | Phase retrieval with signal bias[END_REF]. With the highest SNR convergence is reached, however at a three times higher number of iterations, as can be observed comparing the orange and black curves in Fig. 2(d1), and at a worse final reconstruction (Fig. 2(d2)). The reconstruction without spectral initialization can be improved and accelerated with an exponentially decaying step size (see the green curve in Fig. 2(d1) for the highest SNR case and row (e) for all the noise levels). However, its quality does not reach that of the reconstructions benefiting from the spectral initialization. Notably, an error reduction by a factor of 2 can be observed at the lowest SNR.

Although conceived in combination with GD optimizations, a spectral estimate may be used as initial guess for any iterative ptychographic algorithm. Rows (f) and (g) show the corresponding reconstructions after a PIE algorithm initialized with the spectral estimate and with a flat distribution, respectively. Better reconstructions are obtained with PIE than with GD in the case of highest SNR, as already demonstrated in [START_REF] Hawkes | Springer Handbook of Microscopy[END_REF] and owing to the spatially varying step size [START_REF] Guizar-Sicairos | Phase retrieval with transverse translation diversity: a nonlinear optimization approach[END_REF]. The advantage brought by the spectral estimate is clear in the two cases with lowest SNR, whereas with the cleanest data PIE alone appears to just win out, judging from the reconstruction errors (although the visual difference is more subtle).

All the results shown above were obtained after reconstructing the probe a(x) (plotted in the last row of Fig. 2 with a two-dimensional colormap encoding its amplitude in the lightness and its phase in the hue) from the highest SNR dataset with a regularized PIE (rPIE) algorithm [START_REF] Maiden | Further improvements to the ptychographical iterative engine[END_REF]. Nevertheless, spectral initializations may be implemented in any self-calibrating algorithm which jointly optimizes for both the object and the acquisition parameters. When performing algorithmic self-calibration with alternating projections, the unknowns are optimized one at a time, while leaving the others constant [START_REF] Chen | Quantitative differential phase contrast (dpc) microscopy with computational aberration correction[END_REF]. In this perspective, spectral methods would naturally have a place as the object update step at the earliest iterations of self-calibrating large-scale algorithms using alternating 3)). In (d1) the reconstruction curves for three algorithms at the highest SNR are plotted, while (d2) shows the reconstruction from WF GD. A normalized error metric [START_REF] Maiden | An improved ptychographical phase retrieval algorithm for diffractive imaging[END_REF] is given above each reconstruction, where all reconstructions are compared with a reference reconstruction from an rPIE algorithm, capable of retrieving both the object and the probe [START_REF] Maiden | Further improvements to the ptychographical iterative engine[END_REF]. SI: spectral initialization; GD: gradient descent; WF GD: Wirtinger flow GD; PIE: ptychographic iterative engine. Bottom row, left to right: probe reconstruction, its two-dimensional colorbar, and the phase colorbar used for the reconstructions of the object. Scale bars: 2 mm. projections. They would efficiently drive the solution beyond local minima induced by experimental noise. This procedure would then be followed by the refinement of the acquisition parameters.

(b1) (b2) (b3) (c1) (c2) (c3) (e1) (e2) (e3) (e1) (e2) (e3) (f1) (f2) (f3) (g1) (g2) (g3) (d2) (d1) i WF GD

FIG. 2. Reconstructions of the phase object in

The matrix multiplication involved in the power method implies that the solution is updated after all the ptychograms have been used once. This makes our spectral method a 'batch' algorithm, like gradient-descentbased pytchographic solvers [START_REF] Thibault | Maximum-likelihood refinement for coherent diffractive imaging[END_REF]. On the other hand all the variants of the PIE [START_REF] Maiden | Further improvements to the ptychographical iterative engine[END_REF] are 'serial' algorithms, delivering a new estimate of the solution every time a ptychogram is used. While the former display a better resilience to noise and enjoy parallel computing, the latter typically converge faster [START_REF] Hawkes | Springer Handbook of Microscopy[END_REF][START_REF] Guizar-Sicairos | Phase retrieval with transverse translation diversity: a nonlinear optimization approach[END_REF][START_REF] Yeh | Experimental robustness of Fourier ptychography phase retrieval algorithms[END_REF]. In order to combine the benefits of both classes, a 'stochastic' method based on 'mini-batches' of the full set of ptychograms has been put forward [START_REF] Odstrčil | Iterative least-squares solver for generalized maximum-likelihood ptychography[END_REF] and shown to be particularly well suited to the reconstruction of low spatial frequency phases [START_REF] Odstrčil | Towards optimized illumination for highresolution ptychography[END_REF]. We hereby point out that the power method is fully compatible with this procedure, as already envisioned in the seminal paper [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF].

Finally, note that spectral methods are essentially different from commonly used initialization procedures for iterative ptychographic algorithms. Those are based on synthesizing a low-resolution object from a subset of the measurement data with additional dedicated procedures [START_REF] Tian | Multiplexed coded illumination for fourier ptychography with an led array microscope[END_REF]. Thanks to the similarity between GD iterations and power iterations, we could instead replace early GD iterations with the spectral method, adding no computational burden to the reconstruction framework while, at the same time, leveraging the full dataset to compute the spectral estimate.

In conclusion, we have accelerated and improved ptychographic reconstructions by spectral methods on experimental data. Our results confirm that their application can go beyond random settings, where they were originally developed. Compared to a reference GD algorithm using no spectral methods, three times faster reconstructions were obtained. At different levels of SNR, reconstructions employing spectral methods outperformed standard reconstructions, reaching a 2-fold reduction of the reconstruction error at low SNR. Although our proof-of-principle reconstructions used the simplest form of optimization algorithms for ptychography, we have envisioned the implementation of spectral methods in large-scale optimization frameworks too. This would come at no additional computational cost compared to a GD algorithm. where ā(x -x (l) ) denotes the complex conjugate of a(xx (l) ). Stitching these partial estimates yields the estimate ψ t+1 (x) at the (t + 1)-th iteration.

Appendix D: Different preprocessing functions

Over the years, spectral methods based on different preprocessing functions have been proposed, e.g. T 1 in [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF], T 2 in [START_REF] Marchesini | Alternating projection, ptychographic imaging and phase synchronization[END_REF], and the optimal preprocessing function in the noiseless setting T * in [START_REF] Luo | Optimal spectral initialization for signal recovery with applications to phase retrieval[END_REF], defined as:

T 1 (y) ≡ y T 2 (y) ≡ 0 if y ≤ T 1 if y > T T * (y) ≡ 1- 1 y , ( D1
) where T is a predefined threshold usually set to a quantile of the intensity distribution (for example keeping the top 20% values). The optimal preprocessing function T * was used to obtain the results presented in the main text. Since T * (y) = 1 -y -1 is not bounded below, some eigenvalues of Z may be negative which is detrimental for power iterations. Although we aim at finding the eigenvector associated to the largest eigenvalue of Z, the power method selects the eigenvector corresponding to the largest eigenvalue in modulus. Therefore, care must be taken to avoid large negative eigenvalues. To prevent this, in our reconstructions we set a lower bound -β to the preprocessed intensities. To study in more detail the different spectral methods, we performed a numerical experiment with a Siemens star amplitude-only object (Fig. 3(a)), using an angular spectrum propagator at a Fresnel number of 3.9, a ptychographic overlap ratio of 77% and starting the power method from a random initial guess. Fig. 3 shows different initialization estimates for various choices of the preprocessing function T , the number of power iterations M , and the spatial distribution of the illumination function a(x). We have seen empirically that these three parameters are impacting the most the spectral method performance. Note that when the illumination function is a weakly structured Airy disk (Figs. 3(a-f)), the hypotheses of random sampling required by the spectral method [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF] break down. For this reason, even in a noiseless setting like the one presented in Fig. 3, we cannot expect to indefinitely approach the solution using power iterations, as confirmed by Figs. 3(b,c). Moreover, the optimality of T * is not established with our non i.i.d. random matrix S. However, the preprocessing functions T 2 and T * make the estimates more robust to the power iterations, while delivering informative estimates already after 10 power iterations.

In Figs. 3(h,i), we simulated a random beam with a speckle grain size about 10 times smaller than the illumination shifts, so to boost the diversity of the ptychograms upon translation. This makes the acquisition closer to that of coded diffraction imaging, for which spectral methods were originally developed and have already been applied [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF][START_REF] Mondelli | Fundamental limits of weak recovery with applications to phase retrieval[END_REF]. As a result, a much more informative spectral initialization is obtained, with no need to preprocess the measurement data. Besides an increase of the spatial resolution, we notice a more reliable quantification of the amplitude.

These results can also be seen in the context of ptychography with randomized and structured illuminations, which has been implemented with three main advantages: adding diversity in the ptychograms [START_REF] Odstrčil | Towards optimized illumination for highresolution ptychography[END_REF][START_REF] Stockmar | Near-field ptychography: phase retrieval for inline holography using a structured illumination[END_REF], accessing higher spatial frequencies [START_REF] Bian | Ptychographic modulation engine: a low-cost DIY microscope add-on for coherent super-resolution imaging[END_REF] as in coherent structured illumination microscopy [START_REF] Wicker | Resolving a misconception about structured illumination[END_REF], and reducing dynamic range requirements in X-ray imaging [START_REF] Morrison | X-ray ptychography using randomized zone plates[END_REF][START_REF] Maiden | Soft X-ray spectromicroscopy using ptychography with randomly phased illumination[END_REF].

Appendix E: Amplitude correction for the spectral initialization

Because the spectral initialization is calculated from an eigenvalue problem, one needs to design a strategy to choose its norm, especially for amplitude objects such as the Siemens star (see previous section). Moreover, this normalization is performed pixel by pixel to account for the fact that all the pixels are not sampled uniformly in ptychography.

In order to obtain physically meaningful values of the modulus of the solution, the spectral initialization after M power iterations ψ M was normalized with the spectral initialization ψ ref

Fig. 1

 1 FIG. 2. Reconstructions of the phase object inFig. 1 with different algorithms (rows (a)-(c) and (e)-(g)) at increasing SNR (labels (1)-(3)). In (d1) the reconstruction curves for three algorithms at the highest SNR are plotted, while (d2) shows the reconstruction from WF GD. A normalized error metric[START_REF] Maiden | An improved ptychographical phase retrieval algorithm for diffractive imaging[END_REF] is given above each reconstruction, where all reconstructions are compared with a reference reconstruction from an rPIE algorithm, capable of retrieving both the object and the probe[START_REF] Maiden | Further improvements to the ptychographical iterative engine[END_REF]. SI: spectral initialization; GD: gradient descent; WF GD: Wirtinger flow GD; PIE: ptychographic iterative engine. Bottom row, left to right: probe reconstruction, its two-dimensional colorbar, and the phase colorbar used for the reconstructions of the object. Scale bars: 2 mm.

FUNDING

  Swiss National Science Foundation (200021 160078/1 and P2BEP2 188152); H2020 European Research Council (SMARTIES-724473).

FIG. 3 .

 3 FIG.3. Spectral initializations at different numbers of power iterations, preprocessing and illumination functions (simulated data). Note that, in order for all the images to share the same colorbar, shown on the right, (b) and (c) have been multiplied by 5 and 200, respectively. Scale bar: 1 mm 10 times the wavelength).

  M of a "reference" experiment without object. In other words, ψ ref M is the leading eigenvector ofZ ref = S † diag{T (y ref )}S , y ref = |Sψ ref | 2 , (E1)with ψ ref ≡ [1, . . . , 1] T ∈ R n denoting the transmission function through free space and the superscript T indicating the transpose operation.

ACKNOWLEDGMENTS

LV wishes to thank R. Brönnimann, E. Hack and P. Zolliker (Empa) for fabricating the sample and for contributing constructive comments during the data acquisition.

Appendix A: Intuition of the spectral estimate Spectral methods are simple-to-code algorithms to find initial estimates of the solution to phase retrieval problems, usually combined with gradient descent to refine the spectral estimate. The intuition behind these spectral methods is the following: for one intensity measurement y i = |s † i ψ| 2 (with s † i denoting a row of the matrix S), when s i is correlated with ψ then the measured intensity y i is increased. On the other hand, if s i is orthogonal to ψ, the measured intensity is equal to 0. Therefore we aim to create a matrix dominated by the sampling vectors aligned with ψ. This leads to the definition of a weighted covariance matrix to exploit this information:

which corresponds to Eq. ( 3) in the main manuscript. This matrix is essentially a sum of rank-1 matrices s i s † i , each weighted according to the measured intensity y i . Such matrices first appeared without preprocessing function T , but it was quickly realized how this operation could improve the performance of spectral methods. In the end, T only needs to be an increasing function in order to respect the intuition detailed above. There is a relatively good understanding of spectral methods in the random setting. As the number of measurements p increases, the leading eigenvector of Z converges towards the sought solution ψ. The sample complexity, i.e. how many more measurements p one needs compared to the dimensionality of the problem n, is well characterized when T is bounded above [START_REF] Lu | Phase transitions of spectral initialization for high-dimensional non-convex estimation[END_REF]. In particular, optimal spectral methods actually outperform gradient descent alone, as they require fewer measurements to provide a meaningful estimate [START_REF] Mondelli | Fundamental limits of weak recovery with applications to phase retrieval[END_REF].

Because of this previous heuristic, the recent theoretical advances of spectral methods may be applied to other non-random settings, beyond ptychography. It would be interesting to investigate whether they may apply to other phase retrieval problems, even for example 3D reconstructions.

Appendix B: Comparison with GD update

We would like to emphasize how one power iteration of the spectral method relates to one gradient descent iteration. In particular, as made explicit in the formulas below, the computational complexity of one power iteration is comparable to a gradient descent iteration.

Let us start with the gradient descent update with the amplitude loss function as described in the main text. At each iteration t > 0, the current estimate is computed according to the following equation [START_REF] Yeh | Experimental robustness of Fourier ptychography phase retrieval algorithms[END_REF]:

(B1) On the other hand, with the optimal preprocessing function T * (y) = 1 -y /y, by performing a linear Taylor expansion in y/ y near 1, we obtain the preprocessing function T (y) = 2 y/ y -1 . A power iteration in this case is described by the following map:

Additionally, it is possible to add a term ψ t-1 by replacing Z by Z + I, an operation which does not change the leading eigenvector, but averages the current estimate of the leading eigenvector with the previous one. We thus observe that the two updates are quite similar. Despite the similarities in the equations, note that spectral methods are guaranteed to converge towards the eigenvector of the largest eigenvalue, whereas gradient descent may converge to a local minimum of this non-convex optimization problem.

Appendix C: Non-vectorized version of power iterations

The implementation of the matricial formalism involved after the definition of S typically becomes prohibitively large in a usual ptychographic setting. For example, P is of size Lm × Lm, with Lm the total number of measured pixels, typically much greater than 10 4 . Although this matricial formalism is important to link ptychography with other theoretical settings, it would be unpractical to store P in memory and repeatedly compute multiplications by P to retrieve the leading eigenvector of the associated matrix Z.

Instead, we perform each power iteration by observing that a multiplication by P A corresponds to applying the ptychographic forward model, a multiplication by diag{T (y)} corresponds to an element-wise multiplication with the measured intensities, and applying A † P -1 is equivalent to a backward pass in the ptychographic model (i.e. the backpropagation step of gradient-descent methods). In practice, this operation can be done sequentially, for each ptychogram indexed by l = 1, . . . , L.

Each step provides us with a partial estimate of the solution ψ (l) t+1 (x), which is restricted to the region illuminated at the scan position x (l) . Using the symbols employed for the continuous-domain model in the main manuscript, we obtain:

t+1 (x) = ā(x-x (l) ) P -d T (y (l) (x)) P d {a(x-x (l) )ψ t (x)} , (C1)


domain: spi phys chim math info

text:
P Malbéqui

P.-E Normand

L Budin

G Carrier

SENSITIVITY ANALYSIS OF THE PROPAGATION OF LOW SONIC BOOM TO THE AEROLOGY

Several concepts of supersonic commercial aircraft are emerging and most of them aim at flying over land. The main obstacle for the development of a supersonic aircraft remains the loud and sudden sonic boom felt by the population overflown during the cruise. As low boom designs gain in maturity an acceptable sound level for the sonic booms could be reached. However the acceptability level has still be defined by the regulation authorities as well as the means to evaluate it. The sonic boom sound level on the ground depends on the nearfield signature of the aircraft together with the state of the atmosphere it goes through. This paper examines the influence of the aerology on a low boom signature. The propagation is predicted using the eikonal equation. The aerology is provided by the Integrated Global Radiosonde Archive database, from which one thousand atmospheric profiles altogether have been extracted at 24 locations all over the world. This database contains a representative set of atmospheres, allowing a code-to-code benchmarking. The comparisons done in a collaborative effort with Dassault Aviation and ONERA are in a quite good agreement. Disparities can be explained through the analysis of numerical method implemented in the codes. The sensitivity analysis of the propagation of the sonic boom shows that, despite dispersions on the aerology, the maximum overpressure and the rise time of the pressure signal are close together. In contrast, the extent of the carpet, derived from the cut-off rays on the ground, exhibits a larger deviation.

INTRODUCTION

Unlike procedures used for subsonic aircraft, where ground noise measurements could be backwardpropagated to the aircraft source and then forwardpropagated for prediction in prescribed reference conditions, the sonic boom generated by a supersonic aircraft does not enable us such an approach due to nonlinearities and strong atmosphere effects during long range sound propagation.

As specified in the European RUMBLE project [1], future sonic boom regulation and certification scheme will rely heavily on simulation. As a result, in addition to measurements, the sonic boom prediction will have to play a key role in the future supersonic certification process. Henceforth, prediction capabilities with high level of accuracy is requested regarding both near-field air flow field and ground far-field pressure signature. The RUMBLE activities on the sonic boom prediction capabilities include several items: • to validate the numerical methods for sonic boom prediction, to accurately characterize the near-field signature and predict the sonic boom effects on ground using experimental data from wind tunnel near-field measurements, flight tests data, and codeto-code verification, • to validate wind tunnel techniques for high accuracy near-field signature measurement in cruise conditions, • to develop and validate models for indoor sonic boom prediction accounting for building response, • to design the aerodynamic shape of a low boom demonstrator for the preliminary design of a low boom demonstrator, • to provide recommendations on sonic boom prediction process complying with certification procedures. This paper deals with the prediction of the propagation in the atmosphere from the aircraft to the ground, generated by a low boom concept. The main physical phenomena are addressed using the ray method. The influence of the turbulence in the Planetary Boundary Layer is neglected.

Section 2 summarizes the theoretical background of the ray method, focusing on the techniques and taking into account the specificities of the sonic boom. Comparisons between the codes BANGV and DABANG used by ONERA and Dassault Aviation, respectively, are presented in Section 3, including the configurations of the ISO-Standard atmosphere and the IGRA database. Section 4 describes the sensitivity analysis of the propagation of the sonic boom on the ground, performed through the various conditions available in the IGRA database.

THEORETICAL BACKGROUND ON PROPAGATION OF THE SONIC BOOM

The main specificities of the two codes DABANG and BANGV for predicting the signature are listed below, more details can be found in [START_REF] Dagrau | Simulation de la propagation du bang sonique: de la CFD à l'acoustique non linéaire[END_REF] to [START_REF]BANGV: the Sonic Boom Software, User's Guide[END_REF]. The codes allow for sonic boom propagation under the following configurations:

• Any shape of the aircraft is possible; the near-field signature can be modeled by the so-called Whitham F function or with a cylindrical pressure distribution at a given distance of the aircraft, derived from the CFD or the wind tunnel measurements. • Various aircraft trajectories can be handled, such as cruise, maneuver and acceleration. • A stratified meteorology without turbulence (no temporal fluctuations in the media during the propagation) is assumed, with profiles of temperature, wind, relative humidity and density. The interpolation is done using cubic splines. • A flat and absorbing ground is considered, using the 4 parameters model proposed by Attenborough [START_REF] Attenborough | Acoustical impedance models for outdoor ground surfaces[END_REF].

The two codes are based on the same ray theory, stating that the sound propagating along rays can be determined according to the Fermat principle. This theory is the well-known high frequency approximation, assuming locally plane wave propagation. The amplitude of the wave form is predicted along rays, through the cross-section of the ray-tube, using the Blokhintzev principle. The numerical implementation follows the formulation derived by Candel [START_REF] Candel | Numerical solution of conservation equations arising in linear wave theory: application to aeroacoustics[END_REF]. It is based on the integration of a system of 13 Partial Differential Equations: 4 PDEs to determine the rays, 8 PDEs to estimate the ray-tube cross-section, and 1 PDE for the nonlinear age variable.

To account for the non-linearity of the propagation in the atmosphere, due to the significant magnitude of the peak overpressure in the sonic boom, the codes solve the non-linear Burgers's equation, governing the distortion of the signal, associated with an appropriate technique to automatically handle the shock position. The acoustical absorption due to the thermo-viscous effects and the molecular relaxation of nitrogen and oxygen molecules is performed using a split-step algorithm.

CODE-TO-CODE COMPARISON

Dassault Aviation and ONERA use the codes presented in the previous Section, to validate the capability of the ray method to account for the main physical phenomena affecting the atmospheric propagation of a sonic boom signal, representative of a low-boom aircraft signature. The signal used as the common input is the near-field pressure signal from the NASA C25D low boom concept (see Figure 1). This concept corresponds to the NASA low boom flight demonstrator defined in the 2 nd AIAA Sonic Boom Prediction Workshop. It is the result of the prediction using a Reynolds-averaged Navier-Stokes (RANS) computation at flight unit Reynolds number of 5.7 million per meter. Details of the comparison exercise done at Computational Aero-sciences Branch, NASA Langley and Ames Research Centers can be found in [START_REF] Park | Nearfield Summary and Statistical analysis of the second AIAA Sonic Boom prediction Workshop[END_REF].

The aircraft trajectory is defined as follows: the cruise Mach number is assumed to be constant (no acceleration) and equals 1.6; 0° and 180° heading angles are considered, and the altitude of the aircraft is 15760 meter. The capability to predict the sonic boom primary carpet, as well as the sonic boom signatures over this carpet, is presented for various atmospheric conditions and a given aircraft trajectory, according to the test cases defined in RUMBLE [START_REF] Normand | Far field modeling validation -test cases specifications[END_REF].

The ISO-Standard atmosphere (ICAO 7488/3, 1993) with a constant 70 % relative humidity is first considered. Figure 3 illustrates the temperature and the density profiles. Figure 4 shows a quite good agreement between the two codes of the signature predicted undertrack, both in terms of magnitudes of the positive and negative overpressure of the shocks as well as the rise time. The objective to compare: the cut-off angles, the positive and negative overpressures on the ground, and the metrics: ASEL, BSEL, CSEL, and PL values. The comparison is done by looking at histograms of the absolute difference between the values produced by the two propagation codes. The histogram's number of equally distributed bins is set with the hypothesis that the probability distribution function associated with code-tocode comparison should be Gaussian-like; therefore, the Scott's normal reference [START_REF] Scott | On optimal and data-based histograms[END_REF] rule can be used. For pressure and sound metrics we set the standard deviation to be a third of the expected precision, for cut-off angle comparison this criterion is relaxed because of the scarcity of the results.

Values of the bin-width and expected precision are presented in Table 1. Table 1. Bin-width and expected precision on the predictions by the ray method.

Figure 6 and Figure 7 show the statistical comparisons between DABANG and BANGV for 794 atmospheres, on the cut-off angles, the overpressures, and the metrics, for an aircraft heading east. The expected precisions from Table 1 are represented with vertical red lines. By convention, positive values imply that the BANGV result is greater than the DABANG result, and for negative overpressure and minimum cut-off angle the opposite is implied. Concerning differences in cut-off angle, as much for west heading than east, more than 60% of the time the comparison is within the expected precision. However a non-negligible number of comparison points are largely off precision. Moreover BANGV statistically predicts a narrower carpet. The source of this discrepancy has been identified and comes from the difference in methodology used to determine cut-off angles in the two codes. In DABANG a ray is considered out of the cut-off zone as soon as the ray is refracted upwards (z(t+dt)> z(t), with z(t) the altitude of the ray at time t) and a bisection method is used to determine the maximum and minimum angle for which the ray is not refracted upwards. In BANGV an analytical solution of the cut-off angle is implemented for which in some cases the solution differs largely from the geometrical definition of the cut-off angle. BANGV now implements the two methodologies.

In terms of pressure differences most of the points are outside the expected precision, this is particularly noticeable in the histogram in logarithmic scale shown in figure 8. BANGV predicts most of the time higher positive overpressure and lower negative overpressure. It is interesting to see that this statistical bias is reflected differently amongst metrics. Because of the weighting of the low frequency content from ASEL, BSEL and PL the bias is lowered and the difference distribution tends to be more centered unlike CSEL which maintains the bias towards higher values for the BANGV code. The difference observed is suspected to come from different implementation of the acoustical absorption model. Further investigations are ongoing to determine the difference of absorption and dispersion values given by the model during the ray propagation between the two codes.

We now focus on the cut-off angles for the two locations ANCHORAGE and PRETORIA, where the largest differences between the two codes were observed. While the overpressure differences between DABANG and BANGV do not exceed one Pa, as shown on Table 3 and Table 3, the deviation in the cut-off angles can reach up to 5 degrees. It is here assumed that for these locations, where the impact of the lateral ray on the ground is sensitive to the conditions of the atmosphere, the two different techniques implemented in the codes may explain the differences. In DABANG the impact point is processed using a bisection algorithm, while in BANGV, under downward refracting atmosphere, the emission angle of an horizontal ray can be determined numerically by relating the emission angle, here straightforwardly computed, to the altitude at which the ray remains horizontal [START_REF]BANGV: the Sonic Boom Software, User's Guide[END_REF]

INFLUENCE OF THE METEOROLGY ON THE LOW BOOM SIGNATURE

Previous studies have already analyzed in details the influence of the meteorology on sonic boom propagation.

Along the ground track of a Paris to New York trip, close to those used for Concorde, including an aircraft acceleration [START_REF] Blumrich | Variability of focused sonic booms from accelerating supersonic aircraft in consideration of meteoroligcal effects[END_REF]. Using the Whitham's F-functions, where the source term is adapted from an Airbus mockup of a planned European Supersonic Commercial Transport designed for carrying 250 passengers at Mach 2 in cruise [START_REF] Blumrich | Meteorologically induced variability of sonic boom characteristics of supersonic aircraft in cruising flight[END_REF]; for a hypersonic Mach 6 cruise aircraft [START_REF] Loubeau | Effects of meteorological variability on sonic boom propagation from hypersonic Aircraft[END_REF]. In this paper, the sensitivity of the aerology is dedicated to the low boom concept. It is first detailed for the influence of the aerology at Stuttgart in Germany, where an amount of 56 atmospheres are available, allowing a relevant statistical study. Figure 11 plots the profiles of the magnitude of the wind speed, the temperature, the relative humidity and the density. The relative humidity points out a great day to day variability on the whole runs. The wind speed profiles also show strong variability, with an increase close to 10000 meters, coherent with prevailing westerlies occurring at latitudes between 30 and 60 degrees. In contrast, as expected, the density and the temperature profiles decrease with the altitude with only weak differences between the runs. The temperature decreases typically by 5 degrees °K per km up 12 000 meter, then temperature when reaching the tropopause is roughly constant, and equals 220 °K, due to the disappearance of the atmosphere and the ozone layer. Figure 12 shows the ground carpet pattern (blue curve) and the width of the carpet (red line) of the sonic boom, determined by the two cut-off angles. Under ISO-Standard atmosphere, without speed, the carpet is symmetrical to the trajectory axis of the aircraft, while in the real atmosphere, in presence of a transverse wind, an asymmetry may occur. Then, differences may also appear depending on heading of the aircraft. In general, for each location, the same tendencies as in Stuttgart can be observed: the standard deviation of overpressure and rise time remains small. Moreover comparing the whole locations, except at Seychelles, weak variations of the mean value arise: the overpressure varies from 17 to 20 Pa and the rise time varies in the range 1.4 to 1.6 10 -2 s. In contrast, the variability of the width carpet, which strongly depends on the wind speed, is more significant. While this study is restricted to the aerology encountered during August, its general tendency is in agreement with the one observed in [START_REF] Blumrich | Meteorologically induced variability of sonic boom characteristics of supersonic aircraft in cruising flight[END_REF], where it was observed that the variability is low undertrack, and is higher for lateral angles in terms of carpet width and boom amplitude.

CONCLUSION

A satisfactory agreement on the shape of the signature under track is found, especially the positive and negative overpressure and the rise time, between the two codes DABANG at BANGV, both based on the ray method. Differences in the width of the carpet can be observed for a few locations. This can be explained by the presence of nearly horizontal rays along the ground, where the prediction of the impact is sensitive to the algorithm implemented in each code to find the cut-off rays. This aspect, encountered for a set of singular configurations is investigated by comparing the ray trajectories and the age function. It is assumed that, when solving the same equations, the same results are expected. However for this kind of situations, the assumption of the ray method itself is no more valid. CAA methods based on the Euler's equations taking into account the scattering effect and turbulence of the Planetary Boundary Layer should provide more realistic results, with in contrast, a significant larger CPU time required.

The prediction on the propagation of the sonic boom shows that, despite disparities on the aerology and the climate related to the location, the over pressure and the rise time of the sonic bang on the ground remains close together, with typical standard deviations of a few Pa and 10 -3 s, respectively.

Further analysis using various trajectories including manoeuvers and acceleration should be examined to extend the validity domain of the ray results.

Figure 1 .

 1 Figure 1. Nearfield signature on the aircraft, NASA C25 low boom Concept [8].The nearfield signature used to initialize the ray codes under track (azimuthal angle of zero degree), is plotted in Figure2. The propagation starts at the reduced range R/L = 3, where R is the 330 feet and L is the 110 feet body length.

Figure 2 .

 2 Figure 2. Near field signature under track, R/L = 3, NASA C25 low boom Concept.

Figure 3 .

 3 Figure 3. ISO-Standard atmosphere: temperature and density profiles.

Figure 4 .

 4 Figure 4. Signature on the ground undertrack for ISO-standard atmosphere. Dassault Aviation and ONERA predictions.The atmospheric condition from the Integrated Global Radiosonde Archive available on line is now considered[9]. Data for the whole month of August 2012 from a subset of 24 locations in the world have been assembled by NASA and transformed into a HDF5 Format, to be used by RUMBLE partners. Figure5indicates the 24 locations from which data was extracted from IGRA.

Figure 5 :

 5 Figure 5: 24 Locations all over the word from the IGRA data base.

Cut

Figure 6 .Figure 7 .Figure 8 .Figure 9 .

 6789 Figure 6. Histograms of the difference between BANGV and DABANG, aircraft heading east. Cutoff angles and overpressures.

Figure 10 .

 10 Figure 10. Weighting for metric ASEL, BSEL and CSEL (credit Peter J Skirrow).

Figure 11 .

 11 Figure 11. Set of 56 aerology profiles at Stuttgart, Germany.

Figure 12 .Table 4 .

 124 Figure 12. Definition of the width carpet. Left: ISO-Atmosphere; right: real atmosphere. Table 4 summarizes the prediction of the mean value, ∑ , and the standard deviation, ∑ | | , where xk represents a sample of the positive overpressure, the Rise Time and the carpet width over the N samples at Stuttgart, and for the 0° heading angle of the aircraft trajectory. It is observed that the standard deviation of overpressure and rise time remains weak, typically a few percent of the mean value, while for the width of the carpet it can reach 10 %. over pressure (Pa) Rise Time (s) Carpet width (km) 19.5 0.5 1.4 10 -2 10 -3 55 8Table 4. Mean and standard deviation of the overpressure, rise time and carpet width at Stuttgart.

Figure 13 Figure 13 .

 1313 Figure13finally summarizes the mean values and the standard deviation of the overpressure and the rise time propagated under track, and the width carpet, for the whole locations. In practice the number of samples for a given location typically varies from 10 to 60, so that the total number of samples is about one thousand. For Shiraz, only one sample is available so that no standard deviation is computed. On the right hand side of the graphic, STD and 0-STD bars correspond to ISO-Standard atmosphere results with and without atmospheric absorption, respectively.

Table 2 .

 2 . ANCHORAGE: 50 meter altitude, run 1346328000, August the 30 th 2002, at 12:00.

		min cut-off (°)	Max cut-off angle (°)
	BANGV	-48.53	50.35
	DABANG	-52.43	54.92

Table 3 .

 3 PRETORIA: 1523 meter altitude, run 1343815200, August the 1 st 2002, at 10:00 AM.

ACKOWLEDGMENTS

This project has received funding from the European Union's Horizon 2020 research and innovation programme, under grant agreement No 769896.


Jean Paul Nery

Matteo Calandra

Francesco Mauri

Ab-initio energetics of graphite and multilayer graphene: stability of Bernal versus rhombohedral stacking

come

I. INTRODUCTION

Multilayer graphene has in principle many possible stable configurations. If the lower two layers are fixed in configuration AB, N layers of graphene have 2 N -2 local minima, since each additional layer can take two possible positions. Of all these possible stacking sequences (SSs), the two that are usually observed exhibit periodicity. The most common one has AB stacking and is known as Bernalstacked multilayer graphene (BG). The other one is rhombohedral-stacked multilayer graphene (RG), which has ABC stacking, but its occurence is less frequent 1 . Our interest lies in RG, due to the presence of flat bands close to the Fermi level [2][3][4][5] , which makes it a candidate for exciting phenomena like high temperature superconductivity 6,7 , and also charge-density wave or magnetic states [8][9][10] . In fact, superconductivity might have already been observed in interfaces between BG and RG 11,12 .

The region of the Brillouin zone where the bands are flat increases with the number of layers, up to about 8 layers 4,9 . It has been known for a long time, through X-ray diffraction measurements, that shear increases the percentage of RG 13,14 . Furthermore, we showed recently that shear stress can lead to many consecutive layers of RG 15 . All of the samples where long-range rhom-bohedral order has been observed involved shear stress, due to milling 16 , exfoliation 3,17,18 , or steps in copper substrates used in chemical vapor deposition (CVD) in Ref. 19. In this last work, RG with a thickness of up to 9 layers and area of up to 50 µm 2 was grown, while previous studies were limited to widths of 100 nm (see references within Ref. 19), which is important for device fabrication. Other approaches that could lead to long-range rhombohedral order (and large areas), but that so far have lead to only a few layers of RG, include curvature 20 , twisting 21 , an electric field 22,23 , and doping 23 . According to Ref. 24, anisotropic strain favors BG, so it should be avoided during transfer to a substrate (like hexagonal boron nitride) or upon deposition of metal contacts. The authors also suggested removing all non-rhombohedral parts prior to transfer (via etching, or cutting with the tip of an atomic force microscope 25 ). This is in line with the observation of Ref. 18 that, when transferring flakes with some rhombohedral order, such order was more likely to be preserved after transfer if the whole flake was ABC-stacked.

Despite the immense interest in graphene systems, the relative energies of different SSs are still not well established in general (neither experimentally nor theoretically), and in particular between BG and RG. Understanding the precise external conditions under which RG is favored is fundamental to better control its production and preservation. To achieve this, it would be useful to first determine the energy of different SSs without any such external factors (such as shear, strain, curvature, twisting, pressure, doping, or electric fields). We only know of 2 experimental works that determine the energy difference between RG and BG. More than 50 years ago, Ref. 14 estimated it by using a calorimeter and the enthalpy of formation of potassium graphite, and obtained 0.61 meV/atom. A more recent experiment 23 used small-angle twisted few layer graphene with different doping levels and electric fields. Through the curvature of stacking domain walls, energy differences between -0.05 and 0.10 meV/atom were obtained. The data was fit with a tight-binding model (for the energy dispersion of ABA and ABC) and used E ABC -E ABA = 0.05 meV/atom. Density functional theory (DFT) calculations obtain similar values (although with both positive and negative signs, as we see in the first section in more detail). A related quantity is the stacking fault (SF) energy, which has been measured and calculated to be 0.14 meV/interface-atom (see Table II). Thus, more recent works suggest that the energy differences between BG and RG are in the 0.01-0.10 meV/atom range, and that the value of Ref. 14 is overestimated.

Here, we study the energy differences between BG and RG in bulk, and the energy distribution of all SSs for different number of layers N . Our calculations use DFT within the local density functional approximation (LDA). First, we calculate the bulk (N → ∞) energy differences, and justify why they are so small. Then, it is shown that the low energy configurations at room temperature are probably BG, RG and mixed BG-RG phases with a particular type of interface which we refer to as soft. In addition, the stability of BG and RG is studied at different electronic temperatures and for different amount of layers, and a temperature vs. N phase diagram is obtained. Finally, we introduce a spin Ising model to fit the DFT energies, which helps to better understand how are SS distributed, and can also be used to predict the order of configurations for larger N .

II. BULK ENERGY DIFFERENCES

A. Different functionals and temperatures

We first focus on the bulk energy difference between RG and BG, since there are contradicting results in the literature. DFT calculations have obtained E RG -E BG ∼ 0.1 meV/atom [26][27][28] , but results are probably not well converged (for example, due to the use of a large smearing factor of 0.1 eV [28]). Ref. 29 instead used very dense grids, and obtained E RG -E BG < 0, which agrees with our room temperature result. The electronic temperature can be incorporated by using Fermi-Dirac (F-D) smearing, since it corresponds to the electronic occupation number. Although it is used frequently 9,30,31 , it was not considered in the previous works that calculated E RG -E BG . The result is independent of the smearing function, so different approaches should agree as the smearing parameter δ goes to 0. However, while a finite δ in F-D smearing can be interpreted as the electronic temperature, in other smearing schemes it does not have a clear physical interpretation.

We considered temperatures (δ) of 284 K (which we will refer to as room temperature), 1010 K and 3158 K (0.0018 Ry, 0.0064 Ry, and 0.02 Ry, respectively) and also several functionals (Table I). As we can see, the difference is negative at room temperature (favoring RG), and positive at higher temperatures (favoring BG) for all functionals. Although these energy differences are small, the fact that all results are similar, and that the temperature trend is the same for all functionals, indicates that the results are meaningful. In Fig. 2, we also plot the energy difference as a function of temperature, which shows BG becomes more stable at about 730 K.

In the rest of this work we will use an LDA functional, since it gives good energy differences. For example, the LDA shear frequency ω = 42 cm -1 , related to the curvature of the energy landscape around the minimum AB, agrees well with experiment (see Ref. 15 for more details). Also, Table II shows that it is in good agreement with more sophisticated methods like the adiabaticconnection fluctuation-dissipation theorem within the random phase approximation (RPA) 32,33 and Quantum Monte Carlo (QMC) 34 . The intermediate point between the minima AB and AC corresponds to a saddle point, and LDA agrees well with RPA. For AA stacking, LDA falls in between RPA and QMC. These two last values differ by 41%, so more theoretical studies are needed to precisely determine energy differences in multilayer graphene systems. 35 , an empirical 1/R 6 ij correction is added to the interaction between atoms, where Rij is the distance between them. The results are very similar to the LDA ones, which means that van der Waals corrections are not relevant in this semiempirical form. This is consistent with Ref. 29, where the Grimme correction is of only 0.5 µeV. In the GGA-PBE case 36 , RG is even more favored at room temperature. All functionals give similar results at different temperatures.

The stacking fault energy is of the order of 0.1 meV/interface-atom. LDA and RPA give similar values, but LDA is lower by 0.07 meV/interfaceatom. On the experimental side, Ref. 30 uses anisotropic elasticity theory, revised values for the elastic constants and the experimental data of Refs. 37 and 38, to obtain an average stacking fault energy of 0.14 meV/interface-atom. However, it is an average of 0.05, 0.16, and 0.21 meV/interfaceatom, so there is considerable deviation from the mean value. Thus, further experimental work is needed as well.

B. Magnitude of the differences

The LDA energy difference -0.025 meV/atom at 284 K in Table I corresponds to 0.3 K/atom in temperature units. This suggests that both phases should be observed in the same amount at room temperature, according to the Boltzmann distribution. However, due to the strong intralayer interactions, atoms within one layer move in conjunction. For example, to compare the stability of trilayer ABA vs. ABC, the energy per atom should be multiplied by the number of atoms in the third layer (the whole layer has to move to change the posi- tion from A into C, not just one atom). For a flake of 20 nm, this already corresponds to more than 10 4 atoms, and to energy differences above room temperature. In addition, the exponential character of the Boltzmann distribution translates into large differences in the fraction of each phase (see the Supporting Information, Phase coexistence for small flakes, for further details). This is analogous to the presence of permanent ferromagnetism in cubic Fe, Co and Ni, where the energy of the system depends on the direction of magnetization (property referred to as magnetocrystalline anisotropy) by only ∼ 10 -3 meV/atom (0.01 K/atom) 39 . But domains have thousands or millions of atoms, and the energy difference between a domain oriented in different directions can become much larger than room temperature, favoring one direction over another. In the same way, small energy differences per atom in graphene systems can favor either one of the phases (depending on factors like the temperature or the number of layers), as opposed to getting 50% of each. In making these statements about the stability of BG vs. RG, we are assuming that the difference of the contribution of phonons to the free energy is smaller than the electronic one. The small energy differences between BG and RG are reflected on the fact that RG has been ob-served in both natural and artificial samples 1, 13 . In graphite crystallized by arcing, Ref. 1 observed 80% of BG, 14% of RG and 6% of random stacking. Shear stress has been known to increase the percentage of RG to 20-30% 13,14,40 , and as mentioned earlier, it can lead to long-range rhombohedral order 15 . Other factors like temperature, pressure and curvature also affect which structure is more stable. The energy barrier separating one local minimum from another, of about 1.5 meV/atom, is much larger than energy differences between different graphene SSs, which are of the order of 0.01 meV/atom. Thus, different SSs can be similarly stable, and their proportion can depend significantly on the conditions under which a sample is prepared. III) and 12 layers (Table S1), energies were calculated in (i), (ii) and (iii). b Slice of (a) at a = 2.460 Å. We also calculated the energy difference at the theoretical values of a and c d (for which the structure is fully relaxed), and considered the c d vs. a phase diagram around the theoretical value at room temperature, Fig. 1. For both the experimental and theoretical values RG is favored, and for most of the range of the parameters (see also Fig. S3 in the Supporting Information for the GGA-PBE phase diagram, which always favors RG). This figure suggests that thermal expansion may be relevant when considering higher temperatures. However, when using experimental values of a(T ) and c(T ) [42], the energy difference as a function of temperature does not change much (red dots in Fig. 2), and will be neglected in the rest of this work. Fig. 1 also raises the question of which lattice parameters give energy differences closer to experiments. As indicated earlier, Ref. 1 observed 80% of BG, 14% of RG, and 6% of disordered structures. Samples were produced by arcing (as opposed to exfoliation, which would increase the percentage of RG), indicating that RG is a low energy structure (there are much more disordered configurations, but still account for a lower percentage). Furthermore, Ref. 43 uses exfoliation to obtain few layers graphene and looks at the distribution of SSs for N = 4, for which there are 3 possibilities: ABAB (BG), ABCA (RG) and ABAC (disordered). They observe 85% of BG, 15% of RG, and no disordered structure. They argue that these percentages are similar to those of Ref. 1, and that the complicated pattern of stacking domains cannot be ex-plained by mechanical processing (the stacking of some domains might have still been shifted from BG or disordered to RG though). They conclude that the distribution originates from the pristine structure of graphite. With this in mind, let us look at the N = 4 energy differences at room temperature (Table III) with the parameters of (i), (ii), (iii) and (iv) in Fig. 1(a). Relative to the disordered structure, RG is more stable in (i) (experimental values) and (iv) (theoretical), but is similarly stable or less stable in (ii) and (iii), which suggests the energetics of (i) or (iv) are more representative of experiments. In the rest of the calculations, we will use the experimental values (i). , where E can be the energy of ABAB (BG, the reference state), ABCA (RG) or ABAC (disordered). At the lattice parameters (i) and (iv) (see Fig. 1(a)), RG is favored relative to ABAC, at (ii) they have similar energies, and at (iii) ABAC is favored. In both Refs. 1 and 43, RG is favored relative to the disordered structure, suggesting the experimental parameters in (i) or the theoretical ones (iv) better represent the energetics of the SSs. The experimental values will be used in the rest of the calculations.

The differences in

III. FINITE NUMBER OF LAYERS

In order to better understand the energy distribution of SSs, let us start by looking at T = 284K for different amount of layers N (Fig. 3). For each N , the reference state at 0 meV/atom energy is the corresponding AB-stacked sequence. The SS with the lowest energy is placed in the x axis at 0, the highest energy sequence is placed at 1, and the other configurations are placed equidistantly in between. Degenerate structures are included (like ABAC and ABCB, due to up and down symmetry), to correctly capture the energy distribution. For example, for N = 3 there are only two possibilities, ABA (BG) and ABC (RG), and BG has lower energy. Up to N = 5, BG has the lowest energy. At N = 6, there is a configuration with energy slightly lower than 0, which is RG. It is interesting to note that the lowest energy corresponds to BG or RG for all N , and that the energy distribution seems to converge to a universal curve. a b FIG. 3. a, Energy of all stacking sequences at T = 284 K for several number of layers N , relative to the energy of BG at that N (at each N , BG has 0 meV/atom). The lowest energy state is placed at 0 in the x axis, the highest energy state at 1, and all other states at placed equidistantly in between. There are 2 N -2 states for each N . For N = 3, 4, 5, BG is the most stable structure, while for N ≥ 6 it is RG (see Fig. 5 for a plot of ERG -EBG as a function of the number of layers). In addition, the low energy structures are mixed BG-RG systems with soft interfaces. See Fig. 4 for more details. The energy distribution can be fit with an Ising spin model (see Fig. 7). b Low energy region of a, indicating also the states corresponding to BG (squares) and RG (crosses).

A. Spin notation

A given SS can be mapped to a spin chain, putting a variable σ i with value +1 or -1 at each interface: for a layer with a given position (A, B or C), if the upper layer's position is the next one in the periodic sequence ABC, then σ = +1, and -1 otherwise. See Table IV for some examples. This makes visualization easier, and will be used next for low energy structures, but its main purpose is to later fit the DFT energies with an Ising model.

B. Soft and hard interfaces

To characterize the low energy states, let us notice that there are two types of interfaces between BG and RG, that we refer to as soft and hard. A SS with a soft interface has the form ABAB[AB]CABC. Starting from the left, the structure is BG (up to the bracket ]), while starting from the right the structure is RG (up to the bracket [). Thus, the layers in between [ and ] are both part of the BG and RG sequences. The form of a hard interface is ABABAB][CBACBA, so no layers are shared. The energy of the soft interface is 0.07 meV/atom, while that of the hard interface is 0.18 meV/atom, which is higher, as expected (see Table S2, which also shows the energy of other interfaces). The most and least stable states for N = 12, at T = 284 K and T = 1010 K, can be observed in Fig. 4. At room temperature, the low energy states after RG have soft interfaces. The form + ----is also a soft interface, since it corresponds to A[BA]CB. The only exception is + + + + + + + + -+ +, which would correspond to a stacking fault. The least favored states have several consecutive + signs, followed by consecutive -signs, which correspond to hard RG interfaces. Therefore, the most stable structures at room temperature are BG and RG (together with mixed BG-RG structures with soft interfaces), while random structures are disfavored, in accordance with experiments 1 . This is one of the main results of our work. As the number of layers increases, the energy of RG decreases relative to BG. In the bulk limit (N → ∞), the difference is -0.023 meV/bulk-atom, favoring RG. This value corresponds to the slope of the blue line (T = 284 K) in Fig. 5, which shows the difference between BG and RG for different N at 284 K and 1010 K. At 1010 K, BG is more stable (for any number of layers), which is consistent with experiments. For examples, graphite samples where X-ray diffraction showed presence of both BG and RG had to be heated up above 1300 • C [13] or 1400 • C [14] to start transforming to BG. They changed completely to BG at about 2700-3000 • C [14]. It is interesting to note that the difference in surface energy (the y-intercept) is about 0.08 meV/surface-atom (RG has larger surface energy than BG) at T = 284 K, but it essentially goes to 0 (-0.004 meV/surface-atom) at T = 1010 K. FIG. 5. Energy difference between RG and BG for different number of layers at T = 284 K (blue) and T = 1010 K (orange). The straight lines are a fit for N ≥ 12, and their slope give the energy difference for bulk. At room temperature the difference is -0.023 meV/atom, favoring RG, while at T = 1010 K it is +0.021 meV/atom, favoring BG (which coincide, within error, with the values of Table I).

C. Phase diagram

Since bulk RG is more stable at room temperature, and bulk BG at higher temperatures, a phase transition occurs at some temperature in between. More in general, Fig. 6 shows a temperature vs N phase diagram. Our results are consistent with the fact that growth of few layers graphene typically leads to BG. But they seem at odds with BG being the most common structure at room temperature. However, graphite likely forms in the 350-700 • C range (650-1000 K), which favors BG (Table I, Fig. 2), and under high pressures 44 , which also favors BG according to our calculations (Fig. 1). At higher temperatures, BG becomes more stable (RG transforms to BG, as we mentioned in the previous paragraph), and so it is possible that it remains in this local minimum as it cools down (even if RG is the most stable structure at room temperature). In addition, if BG is the most stable SS for a few layers, due to nucleation, further layers will continue stacking in the AB order. To be more precise, if the starting sequence consists of a few layers of AB stacked graphene + -... + -, then + -... + -+ has less energy than + -... + --. Therefore, it is not unreasonable that beyond certain amount of layers, RG becomes the most stable phase. As mentioned earlier, in our calculations BG is favored up to N = 5, and for N = 4, Ref. 43 also obtained that BG is favored. In terms of the stability for N ≥ 6, Ref. 18 has observed up to 27 layers of RG. But the final structure depended strongly on the direction of transfer (due to the induced shear stress). Another factor that could play a role is doping. According to the measurements of Ref. 23, high pdoping favors RG, while n-doping favors BG. Also, as mentioned earlier, Fig. 5 shows that there is a significant amount of energy in the BG-RG interface, of about 0.08 meV/interface-atom, so defects or dislocations might also favor BG. Thus, it is hard to draw definite conclusions about the stability of pristine BG vs. RG under no external factors.

D. Ising fit

To further understand the energy distribution of SSs, using the spin notation for a system of N layers (N -1 spins), we can consider an Ising type of Hamiltonian,

H = A 0 (N -1) + A 1 N -2 i=1 σ i σ i+1 + A 2 N -3 i=1 σ i σ i+2 + A 3 N -4 i=1 σ i σ i+3 + . . . , (1)
A i corresponds to the i + 1 nearest neighbor (n.n.) interaction energy. In particular, A 0 is just the n.n. (1st) interaction energy and is independent of the SS.

For large enough N , when the interaction between the surfaces becomes small, A i should become constant, since an extra layer just adds an additional A i term. This model is used to fit the DFT energies, as shown in Fig. 7 for N = 12 (the results at T = 1010K can be seen in Fig. S6). For both BG and RG to be favored over other structures, |A 2 | |A 1 | and A 2 < 0, so using only A 0 and A 1 (red dots) gives a very bad fit. Using coefficients up to A 3 , the fit is good, and including A 4 does not improve it much. To keep the model simpler and shorter ranged, we use A i , 0 ≤ i ≤ 3. In the Supporting Information (Figs. S4, S7, and S8), we show that the N = 12 fit also works well for lower and larger N . That is, it maintains the relative order of the DFT energies of the SSs. |A1| and A2 < 0, so N coeff = 2 does not work well. N coeff = 4 works significantly better than N coeff = 3, and is similar to N coeff = 5. A3 is also negative, so sequences with several consecutive + or -(with rhombohedral order) are favored.

IV. CONCLUDING REMARKS

We have shown that for bulk graphite, at the experimental values of the lattice parameter a and interlayer distance c d , RG is more stable than BG for different functionals with electrons at room temperature, while at 1000 K and higher temperatures BG is more stable. Although the energy differences are very small, they are similar for all functionals (and in particular, have the same sign), suggesting they are meaningful and not just arbitrary fluctuations. We then argue, as it occurs in magnetic systems, that the energy difference per atom has to be multiplied by the number of atoms in the sample, resulting in a phase preference at a given temperature (as opposed to 50% of each phase, or even random stacking sequences (SSs) as well). For LDA, we also showed that RG is more stable at the theoretical values of a and c d . A phase diagram of c d vs. a shows that this holds around these values as well, unless for example the interlayer distance is contracted about 3%. This suggests low pressures are better to grow RG. In order to determine the most likely value of c d , we looked at the energetics of N = 4 and compared to experimental data, which showed no disordered struc-ture. While the experimental c d favored RG over the disordered SS, lower values of c d make them equally likely or even disfavors RG. In addition, the distribution of graphite obtained through arcing shows that RG should have lower energy than disordered structures.

We also looked at the stability for a finite amount of layers N , and obtained a temperature vs. N phase diagram. At temperatures of about 700 K or higher, BG is more stable for any N . This is consistent with experiments, since RG disappears when heating up samples with both BG and RG at high enough temperatures. At room temperature, BG is more stable for a few layers, and then RG becomes more stable. Although this seems at odds with BG being the most abundant natural phase, graphite is formed at higher temperatures and pressures, so it could get locked at a local minimum when cooling down. Furthermore, since BG is favored for a few layers, and the interface energy is significant, nucleation could lead to further layers of BG.

We then calculated the energy for all SSs for N = 12, and noticed that the low energy states at room temperature correspond to RG, BG and mixed BG-RG systems with soft interfaces (using the experimental parameters). This could have implications for superconductivity, if indeed it has been observed (or might be observed) in BG-RG interfaces. At high temperatures instead, RG is not favored any more, and the low energy states consist of mostly BG. That is, alternating + and -in the spin notation, except for some consecutive sign, ++ or --. Finally, we showed that the energy distribution can be fit with an Ising model, and that the N = 12 fit can be used to determine the stacking order for larger N (without need to calculate the energy of thousands or millions of SSs).

Determining the very small energy differences in multilayer graphene systems is still a challenging problem, and more experimental and theoretical studies are needed. Our work helps to better understand their stability at different electronic temperatures and for different number of layers, and serves as a starting point to add external factors like curvature, doping or an electric field in firstprinciples calculations.

V. METHODS DFT calculations. Calculations were carried out with Quantum Espresso (QE) 45 , using an LDA functional unless otherwise specified. k-grids for vacuum calculations were done with 80 × 80 × 1 k-grids. For the planewaves, a cutoff of 70 Ry was used in Tables I and III, and in Fig. 1. In the other cases, 60 Ry were used. The cutoff for the charge-density was eight times these values. It is worth pointing out that to converge energies themselves (as opposed to their differences) at the 0.001 meV/atom level, much higher energy cutoffs are needed.

When referring to an interface or surface, energy differences correspond to an area, and are divided by 2 (since there are 2 atoms in the primitive cell), and we write them in meV/interface-atom units. Another common unit in the literature is mJ/m 2 . Otherwise, the energy differences are divided by the total amount of atoms in the structures that are being compared, and the units are meV/atom.

  FIG. 1. a Bulk energy difference as a function of c d (interlayer distance) and a at room temperature. At the experimental and theoretical values, RG is more stable. For low enough values of c d , BG has lower energy, so high pressures should favor BG. For 4 layers (see TableIII) and 12 layers (TableS1), energies were calculated in (i), (ii) and (iii). b Slice of (a) at a = 2.460 Å.

  -+ -+ RG ABCABCABC + + + + + + ++ Soft BG-RG interface ABAB[AB]CAB + -+ -+ + ++ Hard BG-RG interface BABAB][CBAC -+ -+ + ---

FIG. 4 .

 4 FIG. 4. a At T = 284 K, the low energy states are RG, BG and mixed BG-RG states with soft interfaces. The number between brackets indicates the number of soft interfaces. SF stands for stacking fault. b At T = 1010 K, the low energy states have mostly alternating +,signs. The high energy states are similar in both cases, with 2 or 3 consecutive + signs followed by 2 or 3signs (hard RG interfaces).

  FIG.6. Temperature vs. Number-of-layers phase diagram of the most stable SS. At high temperatures, BG is the most stable structure, and also at room temperature for a small amount of layers. At room temperature for N ≥ 6, RG is the most stable structure. However, BG is observed more frequently, and in experiments that try to manipulate multilayer graphene like in Ref.18, the structure seems to more easily transform into BG. On the other hand, since BG is the most stable structure for small N , nucleation could lead to BG for larger N . In addition, natural graphite is formed at higher temperatures, at which BG is more stable, and could remain locked at a local minimum at room temperature. So it remains unclear if RG is the most stable structure for many consecutive layers. The red curve corresponds to a 1/N fit for N ≥ 12 (see the Supporting Information for details).

TABLE II .

 II Comparison between LDA energies (using a = 2.461 Å and T = 284 K for the electronic temperature) and previous works. Energies are relative to BG (ABABABAB for the stacking fault, and AB for the saddle point and AA stacking). In configuration AA, the LDA result (c

d = 3.603 Å) is in between the RPA 33 (a = 2.46 Å, c d = 3.62 Å) and QMC 34 (a = 2.460 Å, c d = 3.495 Å) values. RPA and QMC differ by about 40%, so more theoretical work is needed to determine energy differences more precisely. At the saddle point, LDA (c d = 3.365 Å) agrees well with RPA (c d = 3.42 Å). The stacking fault value of LDA (c d = 3.333 Å) is lower than RPA by 0.07 meV/interface-atom, or about 0.01 meV/atom.

  Table I were calculated using the experimental parameters a = 2.461 Å (inplane) and c d = 3.347 Å (interlayer distance) 41 .

	E RG E BG (meV/atom)	Temperature (K) 50 200 350 500 650 800 950 1100 0.03 0.02 0.01 0.00 0.01 0.02 a and c at T=284 K With thermal expansion

FIG.

2

. Bulk energy difference as a function of temperature. RG is favored at at low temperatures up to about 730 K. At higher temperatures, BG is favored. The blue dots correspond to the experimental latttice parameters at room temperature 41 , while the red dots include experimental thermal expansion values

42

. The results are similar, and thermal expansion will be neglected.

TABLE IV .

 IV Illustration of the spin notation for several SSs.

  6. Temperature vs. Number-of-layers phase diagram of the most stable SS. At high temperatures, BG is the most stable structure, and also at room temperature for a small amount of layers. At room temperature for N ≥ 6, RG is the most stable structure.

		T = 284 K T = 1010 K
	A0	0.0327	0.0423
	A1	0.0019	0.0257
	A2 -0.0497	-0.0288
	A3 -0.0164	-0.0062
	TABLE V. Ising fit coefficients for N = 12, at room
	temperature and 1010 K. At low temperatures, |A2|
	|A1| and both BG and RG have low energies. Also, A2
	is about three times larger than A3, and the low energy
	states have soft-interfaces. At 1010 K, |A1| is larger,
	which disfavors RG.
			However, BG is observed more frequently, and in exper-
			iments that try to manipulate multilayer graphene like
			in Ref. 18, the structure seems to more easily transform
			into BG. On the other hand, since BG is the most sta-
			ble structure for small N , nucleation could lead to BG
			for larger N . In addition, natural graphite is formed
			at higher temperatures, at which BG is more stable,
			and could remain locked at a local minimum at room
			temperature. So it remains unclear if RG is the most
			stable structure for many consecutive layers. The red
			curve corresponds to a 1/N fit for N ≥ 12 (see the
			Supporting Information for details).

  FIG. 7. Fit of the DFT energies using the Ising model, Eq. (1), for N = 12 layers at T = 284 K. RG is the lowest energy state, and by definition, EBG = 0. To favor both BG and RG, |A2|

		0.035
	E Ising fit (meV/atom)	0.005 0.000 0.005 0.010 0.015 0.020 0.025 0.030		N coeff = 2 N coeff = 3 N coeff = 4 N coeff = 5
		0.010
		0.015	0.01	0.00 E DFT (meV/atom) 0.01 0.02	0.03

ACKNOWLEDGMENTS

We acknowledge support from the European Unions Horizon 2020 research and innovation programme Graphene Flagship under grant agreement No 881603. M.C. also acknowledges support from Agence nationale de la recherche (Grant No. ANR-19-CE24-0028). The authors would like to thank Stony Brook Research Computing and Cyberinfrastructure, and the Institute for Advanced Computational Science at Stony Brook University for access to the high-performance SeaWulf computing system, which was made possible by a $1.4M National Science Foundation grant (#1531492).