Elastic propagation of fast electron vortices through amorphous materials

This article describes the scattering behavior of electron vortices inside amorphous samples. It focuses on the vortex purity, net angular momentum transfer, and statistical variations due to random beam and atom positions.


Introduction
The study of electron vortex beams (EVBs) is a highly active field of research in the context of transmission electron microscopy (TEM). Pure EVBs are characterized by a phase distribution proportional to the azimuthal angle and can thus be written in the form ðr; 'Þ ¼ f ðrÞ expðim'Þ, where ðr; 'Þ are polar coordinates, f is the radial amplitude function and m is the so-called topological charge. It can directly be verified that EVBs are eigenfunctions of the orbital angular momentum (OAM) operatorL L z withL L z ¼ hm , which implies that EVBs carry angular momentum and, by extension, a magnetic moment of m B (Bliokh et al., 2011).
The fact that these EVBs carry OAM has led to the proposition and demonstration of many applications ranging from the measurement of magnetic properties with atomic resolution (Verbeeck et al., 2010;Rusz et al., 2014;Idrobo et al., 2016;Schachinger et al., 2017), the study of the dynamics of Landau states (Schattschneider, Schchinger et al., 2014;Schachinger et al., 2015), sample chirality (Juchtmans et al., 2015) and symmetry properties of plasmon resonances (Guzzinati et al., 2017), to the manipulation of nanoparticles . Despite the huge potential of EVBs and the fact that their creation and propagation through vacuum are well understood (Schattschneider & Verbeeck, 2011;Schattschneider et al., 2012;Schachinger et al., 2015), knowledge of their propagation through matter is still somewhat lacking. This is especially surprising since earlier studies showed that elastic scattering in crystals can drastically change the OAM of the beam (Lö ffler Xin & Zheng, 2012;Lubk et al., 2013).
Particularly important -and little investigated -is the propagation of EVBs through amorphous materials. Firstly, such materials are used increasingly often for producing EVBs by means of specially designed phase masks (Harvey et al., 2014;Shiloh et al., 2014;Grillo et al., 2014). Secondly, they are a common support, e.g. for nanoparticles. Thirdly, EVBs would allow techniques such as energy-loss magnetic chiral dichroism (EMCD) for measuring magnetic properties down to the nanoscale in crystalline samples to be applied also to amorphous materials (Schachinger et al., 2017). However, it is usually assumed that an as-produced, ideal vortex beam stays that way and propagates practically unperturbed through the sample. Whether or not that is the case and, if so, to what extent is studied in this work.
This paper is structured as follows: first, we give a brief overview of the theory in Section 2. To that end, we rewrite the multislice approach used throughout this work in a cylindrical coordinate system suitable for the analysis of EVBs. From that, we deduce some general statements about the propagation behavior of EVBs. In Section 3, we give a detailed account of the numerical simulations performed in this work. In Section 4, the results of the numerical simulations are presented, which are subsequently discussed in Section 5.

Theory
The starting point for describing the propagation of electrons through matter is Schrö dinger's equation. Throughout this work, we will adopt a paraxial multislice approach (Kirkland, 1998). In this approach, the sample is cut into thin slices and the propagation of an electron wavefunction through slice n is given by where r ? is the 2D coordinate vector in the x-y plane perpendicular to the main propagation direction z, nÀ1 ðr ? Þ is the wavefunction incident on the nth slice, n ðr ? Þ is the wavefunction exiting the nth slice, k z is the z component of the wavevector, t n is the thickness of the nth slice,4 4 is the Laplace operator, is the so-called interaction parameter and v z;n ðr ? Þ is the electrostatic potential of the slice projected along the z direction. In equation (1), the exp½Àiv z;n ðr ? Þ term describes (instantaneous) elastic scattering, while the exp½ðit n =2k z Þ4 4 term describes the free-space Fresnel propagation through the slice. To propagate the electron beam through the entire sample, many such individual propagation steps have to be performed. Note that4 4 and v z;n ðr ? Þ do not generally commute, so the exponentials cannot easily be reordered.
Here, we are primarily interested in the evolution of the different OAM components, so we expand the terms in equation (1) into the eigenstates expðim'Þ of the OAM operatorL L z ¼ Àih -@=@': n ðr ? Þ ¼ P m f n;m ðrÞ expðim'Þ exp½Àiv z;n ðr ? Þ ¼ P V n; ðrÞ expði'Þ; with f n;m ðrÞ ¼ 1 2 where ðr; 'Þ denote the polar components of r ? . In physical terms, m denotes the topological charge of a vortex component with an OAM of mh -. With these definitions, equation (1) reduces to X m f n;m ðrÞ expðim'Þ i.e. the elastic scattering transforms the set of radial components ff nÀ1;m ðrÞg m 7 !fg n;m ðrÞg m . The action of the Laplacian operator, i.e. the Fresnel propagation between the slices, is best viewed in reciprocal space. There, the 2D Laplacian reduces to jk ? j 2 and the OAM distribution is maintained , giving X m f n;m ðrÞ expðim'Þ where ðk; ' k Þ are the polar coordinates of the vector k ? , F k ? !r ? denotes the 2D Fourier transform from reciprocal to real space and J m is the Bessel function of the first kind of order m. It can be seen that the redistribution of intensity between different OAM components happens due to the elastic scattering in the electrostatic potential v z [see equation (3)], while during the Fresnel propagation, only the radial distributions evolve but no intensity is transferred between different OAM components. The potential scattering term can also be written in vector form as g n ðrÞ ¼ V n ðrÞ Á f nÀ1 ðrÞ; where ½V n ðrÞ m;m 0 ¼ V n;mÀm 0 ðrÞ is a Toeplitz matrix, i.e. a matrix in which the values along each diagonal are constant, or, equivalently, where adjacent columns (and rows) are identical apart from a shift by one element.
There are several noteworthy points here. First of all, scattering from a component m to a component m þ m takes place only if there exists some r for which V n;m ðrÞ and f nÀ1;m ðrÞ are both non-negligible. On the one hand, this reflects the obvious fact that only those areas of the potential affect the beam in which the beam intensity is non-vanishing. On the other hand, it also implies certain symmetry properties (see Section 2.1).
Secondly, one can expect m ¼ 0 to be the dominant term for thin slices. This results from the fact that, for thin slices, v z is small. Thus, the potential can be written in weak-phaseobject approximation as exp½Àiv z;n ðr ? Þ ' 1 À iv z;n ðr ? Þ; showing that there is a large constant term, which results in a large m ¼ 0 contribution.

Symmetry constraints
Symmetry plays an important role in the scattering behavior of electron beams, especially in crystalline specimens. Even though the potential typically does not exhibit strict symmetries in amorphous materials, it can still show certain 'approximate' symmetries, i.e. atomic arrangements that deviate only slightly from a symmetric case. In fact, while in crystalline samples symmetries typically only hold for certain special, high-symmetry points such as atomic columns and are severely broken if the electron beam is positioned off-column, the random distribution of atoms in amorphous systems means that the same symmetry properties hold in an approximate sense fairly independently of the beam position. Thus, a closer investigation of the symmetry constraints for OAM transfer seems worthwhile. We want to emphasize that this subsection does not only pertain to amorphous materials but also to crystalline ones.
Here, we consider the inherently 2D case in the plane perpendicular to the beam axis (i.e. in a slice). More precisely, we study the transformation properties of the potential scattering term exp½Àiv z;n ðr; 'Þ under the point group Oð2Þ, which contains rotations and reflections (as well as arbitrary combinations of them).
For the case of rotations, we assume that the potential has a -fold rotational symmetry, i.e. v z;n ðr; ' þ 2=Þ ¼ v z;n ðr; 'Þ. Inserting this into equation (2) yields using the summation formula for finite geometric series. Therefore, in the case of a -fold rotational symmetry of the potential around the beam axis, V n; 0 8 = 2 Z, i.e. intensity can only be redistributed between OAM components which differ by an integer multiple of h -.
For the case of reflections, we assume that the potential is symmetric with respect to a mirror line inclined by an angle ' 0 with respect to the x axis, i.e. v z;n ðr; Since the cosine is a symmetric function, it follows that, in the presence of a reflection, V n; ðrÞ ¼ expðÀ2i' 0 ÞV n;À ðrÞ, i.e. the þ and À components differ only by a phase factor. The case in which the scattering coefficients for þ and À components have the same absolute value may lead to the hypothesis that, in such a case, no net OAM can be transferred as both scattering events happen with the same probability. However, this hypothesis clearly cannot be true as an arbitrary potential exhibiting only a mirror symmetry is not circularly symmetric and hence does not commute with the Hamiltonian. Therefore, Heisenberg's equation of motion together with Ehrenfest's theorem dictate that the net OAM has to change over time. The solution to this conundrum lies in interference effects.
While the train of thought of equal probabilities is correct for single scattering, it breaks down when considering multiple scattering as depicted in Fig. 1 (for ¼ AE1). There, it is clearly visible that, after a single potential scattering event in the first slice, the m À 1 and the m þ 1 components have the same total intensity even though their phase structure is obviously different. The propagation behavior of the two components is also different, owing to the different orders of Bessel functions in equation (4). However, as the Fresnel operator is unitary, the total intensity does not change during propagation.
The situation is different after the second slice, though. After the second potential scattering event, the m þ 1 component is given by the coherent superposition of two contributions. The first one stems from the portion of the beam that was first scattered with m ¼ 1, then propagated as m þ 1, and then scattered again with m ¼ 0. The second one stems from the portion that was first scattered with m ¼ 0, then propagated as m, and then scattered with m ¼ 1. The situation for the m À 1 component is analogous, but not identical. Since propagation and potential scattering do not commute and the propagation is dependent on m, the interference patterns emerging from the coherent superpositions can be different for the m À 1 and the m þ 1 components, thus leading to different intensities of the two components as indicated in Fig. 1. This, in turn, leads to a change of the OAM expectation value and, hence, to a net transfer of OAM, even though each individual potential scattering event is (quasi-)symmetric in amplitude for positive and negative m.

Radial dependence
Another interesting question is how the OAM transfer depends on the radius, which translates into the question of how the behaviors of smaller and larger beams differ. For increasing r, larger and larger OAM transfers will become important. In fact, it is reasonable to assume that the dominant OAM transfer mh -6 ¼ 0 should scale proportionally to r. This can be deduced by comparing the mean atomic distance a with the circumference of a circle with radius r. Since the mean distance is constant throughout the sample but the circumference scales linearly with r, the ratio of the two scales as 1=r. For large r, this can be seen in a very crude approximation as the average period p $ a=ð2rÞ of a periodic oscillation of the potential as a function of '. Thus, the frequency of this oscillation (which corresponds to the OAM transfer) is proportional to 1=p / r. Consequently, one can expect that larger OAM transfers become more important with increasing beam size.
Obviously, for large r, there is also more room for variations, i.e. deviations from a perfect periodic oscillation with period p. Therefore, it can also be expected that the spread of possible OAM transfers should increase with increasing r. As an alternative argument leading to the same conclusion, one can invoke the uncertainty principle ½'½L z $ const: Schematic of the evolution of different OAM components upon transmission through two slices. The second column represents the (intially pure) vortex of order m (the images show m ¼ 1), the first column shows the m À 1 component, the third column shows the m þ 1 component, the fourth column shows the total wavefunction for reference, and the fifth column shows the scattering potential (where applicable). Blue arrows depict potential scattering while red arrows indicate Fresnel propagation (P P). For clarity, propagation distances are exaggerated and simple single-atom model potentials are used. Dashed arrows symbolize additional scattering contributions that are omitted here. The numbers in the first and third columns give the components' intensities (with the total wavefunction normalized to 1). The index for the slice number n and the coordinates r; ' were omitted. Brightness signifies amplitude, color signifies phase as depicted in the inset in the top left. angular extent of the order of ½' $ a=r. Thus, one can expect the standard deviation of the OAM, ½L z , to scale roughly proportionally to r as well.

Expectation value
For some applications such as nanoparticle manipulation, the individual components of the OAM play only a secondary role compared with the expectation value of the OAM operatorL L z , which corresponds to the total, net OAM of the beam. Directly in front of the nth slice, this expectation value is given byL while behind the slice it is given bŷ where the last equality holds due to Parseval's theorem (namely that the integrals over the absolute value squared of a function and its Fourier transform are identical), or, equivalently, due to the closure relationship of Bessel functions. Clearly, hL L z i is affected by asymmetries in the OAM component intensity distribution. From the derivations above, one can therefore expect that hL L z i changes more for larger beams (which produce a wider spread of OAM component intensities) as well as for strongly scattering materials (for which jV n; j is relatively large for 6 ¼ 0).

Numerical simulations
We performed extensive numerical simulations for three amorphous model systems: Si 3 N 4 , which is commonly used as support material and for phase masks; Pt, which is commonly used as a focused ion beam (FIB) protection layer and in absorption masks; and Fe 0.8 B 0.2 , a magnetic material used, e.g., in transformers, which could be interesting for EMCD. All simulations were carried out multiple times (see Table 1) for randomly different atom arrangements to get an idea of the variations of the various results. The simulations were performed in Cartesian coordinates using a multislice code based on the one described by Kirkland (1998). Likewise, the atomic potentials were also taken from Kirkland (1998).
For all samples, a 100 Â 100 Å area was simulated with 512 Â 512 pixels using thicknesses in the range of 0 to 500 Å with a slice thickness of 2 Å . All simulations were performed with a 200 keV incident beam which initially was in an OAM eigenstate with L z ¼ h -. The convergence angles were in the range of 1 to 25 mrad, corresponding to waist radii in the range of approximately 10 Å to 0.4 Å (see Fig. 2). For the sake of straightforward interpretation, the experimental conditions were assumed to be ideal, i.e. the microscope lenses were assumed to be perfectly aberration-corrected and no broadening due to a partially incoherent source or motion of the atoms was included.
The atomic positions were generated at random, taking care that the overlap between adjacent atoms was as small as possible (i.e. rejecting atoms that were too close to already placed atoms). The used densities are summarized in Table 1.
To evaluate the OAM components, the resulting wavefunctions n ðr ? Þ after each slice were first transformed to a polar representation n ðr; 'Þ using a fixed ðr; 'Þ grid with 256 Â 1024 pixels. Then, the transformation '7 !m was carried out by separately Fourier-transforming each line of constant r, yielding n ðr; mÞ. Finally, the result was summed over the radius to obtain the total intensities  Table 1 Densities of the materials and number of random configurations used in the simulations.
The mass densities were used as reference. The atom densities were the ones used in the simulations.

Figure 2
Dependence of the m ¼ 1 beam waist radius r on the convergence semiangle for 200 keV electrons (Lö ffler, 2013). the physical quantities due to the fact that no two samples and no two positions on a sample are identical. Fig. 3 shows some examples of the data produced by the simulations and during the analysis. In particular, it shows that, as predicted, the redistribution of intensity between different OAM components produced by the scattering potential is approximately symmetric but the resulting wavefunction has a distinctly non-symmetric OAM component distribution around the initial m ¼ 1 component. Under the given conditions, the m ¼ 1 component still exhibits the highest intensity I 1;n ' 0:15, but the components m 2 fÀ2; À1; 0; 2; 3g have considerable intensities of the order of I 1;n =2. Therefore, their sum greatly exceeds I 1;n . Even higher orders in the range À20 < $ m < $ 20 also have sizable intensities, and even higher orders (not depicted) still have tiny contributions, further emphasizing the broadness of the m distribution. Interestingly, though not surprisingly, different m components contribute strongly at different radii. In addition, the theoretically predicted linear increase of both the domi-nant m 6 ¼ 0 contributions as well as the m spread in the scattering potential are clearly visible. Fig. 4 shows the dependence of several key quantities on the convergence semi-angle (which is related to the beam size, see Fig. 2) of the incident beam as well as the thickness of the sample for the three simulated systems.

Results
The most striking property is that while the overall features of the graphs are comparable between the three different systems, the numerical values differ greatly. Taking the maximum ½L L z as an example, it changes from '16 for Si 3 N 4 to '36 for Fe 0.8 B 0.2 to '66 for Pt. A similar trend is visible for hL L z i. This phenomenon correlates nicely with the mass density of the three systems. Even though the atom density is comparable for Si 3 N 4 and Fe 0.8 B 0.2 and is lower for Pt, the mass density increases from Si 3 N 4 to Fe 0.8 B 0.2 to Pt (see Table 1), owing to the fact that Pt atoms are much heavier than, e.g., Fe atoms. Since heavier atoms generally scatter more strongly, it is logical that such systems produce stronger OAM deviations.
With respect to the changes of the expectation value hL L z i, Fig. 4 shows that the largest net OAM transfers occur for small and large convergence angles and the smallest deviations occur typically around the range 5-10 mrad, especially for Some examples of the data produced by the simulations and during the analysis. Left: wavefunction at the incident plane of the sample, z ¼ 0; center: wavefunction at the depth of z = 500 Å ; right: scattering potential expðiv z;n Þ of one slice. First row: data in Cartesian coordinates; second row: data in polar coordinates; third row: intensity of the expðim'Þ components as a function of m and r; fourth row: total OAM intensities integrated over r. For the wavefunctions shown in the first two rows, the phase is displayed as color (see Fig. 1) and the amplitude is displayed as brightness. For the scattering potential in the first two rows of the right column, the argument of the complex exponential is shown. The incident beam was a pure m ¼ 1 vortex with a convergence semi-angle = 10 mrad incident on the amorphous Pt sample. In all cases, only a subset of the entire data set is shown and the contrast was enhanced to improve visibility. small to medium thicknesses. This can be related to the size of the beam as it propagates through the sample. For small , already the incident beam is large compared with interatomic distances and it stays that way all throughout the sample. Thus, large m are possible from the very beginning of the propagation. For large , the diameter of the incident beam is small, but the beam size increases considerably during propagation. Thus, although initially only small m are viable, larger and larger m become dominant as the beam propagates further through the sample. Conversely, a beam with a mid-range represents a good compromise between small initial size and small growth during propagation, thereby restricting the maximal significant m and, consequently, the variation of hL L z i. A similar result was also found for classical EMCD (where vortex beams are generated during inelastic scattering and subsequently analyzed interferometrically) in crystalline samples (Lö ffler & Hetaba, 2018).
For ½L L z , i.e. the OAM uncertainty or, equivalently, the spread of the m distribution, the picture is very similar. Small lead to a very large increase in with thickness. For medium in the range of 7 to 13 mrad, is smallest, while it increases again for large .
Interestingly, the dependence is different for the m ¼ 1 intensity, which gives an indication of 'how much' of the original, incident-beam structure actually is present at a given thickness. Fig. 4 shows that I 1;n obviously decreases with thickness, but is mostly independent of . In other words: even though the net OAM and the m distribution depend strongly on the beam size through the convergence angle and although there is complex multiple scattering going on back and forth between different m components at different radii (as visible from Fig. 3), the overall intensity of the m ¼ 1 component seems to be fairly predictable.
To investigate the intensity of different m components as well as the expectation value in more detail, Fig. 5 shows graphs of the intensity of m ¼ 1 as well as the adjacent components m ¼ 0; 2 and hL L z i for different convergence semiangles as a function of thickness. The adjacent components were selected because, for applications that depend on the fact that the beam is in an m ¼ 1 eigenstate (such as, e.g., EMCD), typically close-lying other components are more difficult to separate than far-removed ones. As an example, an m ¼ 100 vortex would have practically zero intensity everywhere where an m ¼ 1 vortex is strong, thus making it easy to separate and block, e.g., by an aperture.
As before, the overall behavior of the curves is roughly similar for the m components of all three studied systems, except for the scale of the thickness dependence, which, again, is more dramatic for heavier specimens. Nevertheless, there are several noteworthy aspects visible in the graphs. In the first several å ngströ ms, the decay of the m ¼ 1 intensity as well as the increase of the adjacent components are practically linear. This is to be expected as, for a dominant m ¼ 1 component, the transitions 17 !0 and 17 !2 will be much more probable than the scattering 07 !1, 07 ! À 1 etc. However, after several å ngströ ms, all depicted components start to deviate from their linear behavior. The m ¼ 1 intensity decrease starts to slow down as soon as it reaches '70% of  the initial intensity, while the m ¼ 0; 2 intensities become almost constant somewhere in the range of 10-30%. As the thickness increases, the m ¼ 0; 2 components seem to asymptotically tend towards a similar behavior as the m ¼ 1 component, as is visible for Pt at = 25 mrad and, to some degree, already at = 10 mrad. Note that, in all cases, the variation over several runs clearly indicates that the results are statistically significant, although the variability naturally is larger for larger mass density.
As already shown in Fig. 4, the decrease in m ¼ 1 intensity does not depend strongly on the convergence angle. However, the increase of the adjacent components is influenced by the convergence angle. At the same time, the statistical uncertainty increases for increasing convergence angles (i.e. smaller beam waists). This can be attributed to the fact that for sufficiently small beams (i.e. smaller than the interatomic distance), the propagation behavior is crucially dependent on the (random) relative position of the beam with respect to close-by atoms, whereas for large beams, the effect is averaged over many atoms.
Another interesting result can be found in the behavior and statistical variation of the expectation value hL L z i. For Si 3 N 4 , the deviation from his marginal and fairly well contained in the statistical error. For heavier systems, the deviation from hbecomes much stronger -with a general trend towards decreasing hL L z i -but also the statistical variation between different simulations increases dramatically.

Discussion
Whether the results presented here are encouraging or discouraging depends on the application at hand, the system under investigation, and the chosen experimental parameters.
If pure vortex beams are required, low mass densities as in the case of Si 3 N 4 and low thicknesses are definitely preferable in order to retain a high intensity in the m component of the incident beam as well as little variation for different atom configurations. This also implies that holographic phase masks fabricated on a thick Si 3 N 4 membrane can be subject to a considerable loss of mode purity.
If a high net OAM transfer is sought (e.g. in the case of nanoparticle manipulation), high mass densities as in the case of Pt as well as thick samples and medium convergence angles should be used. This ensures a large OAM transfer while retaining moderate statistical variations for different atom positions.
It should be noted that real sample densities, interatomic bond lengths and scattering strengths may differ from the ones presented here, e.g. due to the use of different materials. In addition, the sample density is influenced by deposition and Expectation value hL L z i and total intensities of the m ¼ 0; 1; 2 components for Si 3 N 4 (top), Fe 0.8 B 0.2 (middle) and Pt (bottom) for = 1 mrad (left), = 10 mrad (center) and = 25 mrad (right) as a function of thickness z. The intensities of the m ¼ 0 and the m ¼ 2 components have been magnified by a factor of 5 as indicated in the legend to improve visibility. The shaded areas indicate one standard deviation as derived from multiple simulations. preparation parameters. However, the simulations presented in this work span from fairly low to quite high mass densities and scattering strengths, thus giving a general insight into how other samples will behave in general.
All the simulations presented in this work were performed under ideal conditions, including no incoherent source size broadening (ISSB), no atom movement and no lens aberrations. Both ISSB and atom movement would lead to an effectively different relative position between the beam and the atoms for each electron in the beam. This is conceptionally equivalent to the averaging over several random atom configurations as done in this study. Lens aberrations generally lead to a coherent broadening of the beam compared with the ideal case. While in such a situation the details of the amplitude and phase of the beam change, the overall results should be the same as those presented here when considering the appropriate beam size (see Fig. 2).

Conclusion and outlook
In this work, we presented extensive simulations of the propagation of electron vortex beams through amorphous materials. To that end, we have rewritten the multislice approach in cylindrical coordinates to get some theoretical insight into the vortex propagation, such as the beam-size dependence of the redistribution of intensity between different m components and the possibility of net OAM transfer despite the fact that the probabilities for transferring AEmhare (approximately) equal. In addition, we have also described the influence of the point-group symmetry on the vortex propagation.
The numerical simulations were performed for the three amorphous model systems Si 3 N 4 , Fe 0.8 B 0.2 and Pt for a wide range of convergence semi-angles and thicknesses. Besides corroborating the theoretical results, the numerical data allowed us to quantify the net OAM transfer, the spread of vortex components that is related to the uncertainty principle and, thus, the purity of a vortex state, as well as the intensity behavior of the most important vortex components. The results showed that in order to retain high purity upon propagation, low-mass-density samples with small thickness should be chosen, while large net OAM transfers can best be achieved in heavy, thick samples. In both cases, intermediate convergence semi-angles around $ 10 mrad proved beneficial.
The results presented in this work will allow theoreticians and experimentalists alike to choose the material for their studies with electron vortices more efficiently. Although this work does not completely replace full simulations for future studies, it does give some general insight into the propagation behavior of EVBs and makes predictions for a large range of systems and experimental parameters. As such, it promises to contribute to future enhancements not only of the fabrication but also of the applications of EVBs.

Funding information
TS acknowledges financial support by the Austrian Academy of Science (Ö AW) for the DOC scholarship and the 'Hochschuljubilä umsstiftung der Stadt Wien' (project H-294689/ 2016). The authors acknowledge TU Wien University Library for financial support through its Open Access Funding Programme.