Equilibrium structures of anisometric, quadrupolar particles confined to a monolayer

We investigate the structural properties of a two-dimensional system of ellipsoidal particles carrying a linear quadrupole moment in their center. These particles represent a simple model for a variety of uncharged, non-polar conjugated organic molecules. Using optimization tools based on ideas of Evolutionary Algorithms, we first examine the ground state structures as we vary the aspect ratio of the particles and the pressure. Interestingly, we find, besides the intuitively expected T-like configurations, a variety of complex structures, characterized with up to three different particle orientations. In an effort to explore the impact of thermal fluctuations, we perform constant-pressure Molecular Dynamics simulations within a range of rather low temperatures. We observe that ground state structures formed by particles with a large aspect ratio are in particular suited to withstand fluctuations up to rather high temperatures. Our comprehensive investigations allow for a deeper understanding of molecular or colloidal monolayer arrangements under the influence of a typical electrostatic interaction on a coarse-grained level.


I. INTRODUCTION
In recent years there is increasing interest in understanding the film morphologies of anisotropic, conjugated organic molecules at inorganic surfaces. A prime example are systems of optically active molecules at inorganic surfaces; these so-called hybrid inorganic/organic systems (HIOS) represent a very promising material class in optoelectronics. [1][2][3][4] Typically, the corresponding organic molecules are strongly anisotropic in shape and are characterized by complex charge distributions, often dominated by a quadrupole moment. 5 By manipulating the orientational structure of such organic layers it is possible to tune the efficiency of the charge carrier transport 6 and thus to optimize the efficiency of the hybrid system. Somewhat earlier, organic molecules at interfaces have attracted attention in the context of so-called Langmuir monolayers, that is, two-dimensional (2D) films of typically amphiphilic molecules constrained to a liquid-gas interface. 7,8 Again, the orientational and translational structure in such systems can be quite complex.
From the theoretical side, a full microscopic treatment of an organic molecular layer is still challenging, and this holds particularly for HIOS (where the substrate is typically patterned). On the level of atomistic molecular a) Electronic mail: thomasheinemann@mailbox.tu-berlin.de b) Electronic mail: moritz.antlanger@tuwien.ac.at c) Electronic mail: martial.mazars@th.u-psud.fr d) Electronic mail: klapp@physik.tu-berlin.de e) Electronic mail: gerhard.kahl@tuwien.ac.at dynamics (MD) simulations, one has to consider a large number of interactions such as atomic bonding, van der Waals (Lennard-Jones) potentials, contributions due to bending or torsion of a molecule (see e.g. Ref. 9), and electrostatic interactions. Even then, it often turns out that atomistic MD calculations do not correctly predict the arrangement of several particles, one main reason being the insufficient treatment of polarizability effects. 10 An example that demonstrates this issue is the quadrupolar molecule benzene, where a representation via atomic point charges 11 or multipoles 12 fails to reproduce the corresponding dimer configuration. This suggests to employ quantum chemical approaches to evaluate the electronic structure of the entire system, instead. However, at the moment, this level of complexity is computationally far too demanding, especially at finite temperatures.
Given these difficulties, the goal of the present paper is to understand generic aspects of the structural behavior of quadrupolar, anisotropic molecules at interfaces based on a coarse-grained model. Specifically, we investigate a 2D many-particle system composed of ellipse-shaped particles with an embedded, axially symmetric quadrupole tensor (higher-order multipoles are neglected). The quadrupole is oriented along one of the ellipsoidal axes as it is the case, e.g., in benzene 13 or naphtalene, whose quadrupole tensor nearly has axial symmetry. 14 We further assume that the center of the particles is confined to a 2D plane, and that the molecules are allowed to rotate only within this plane.
Interestingly, 2D quadrupolar systems are, so far, rather unexplored. This contrasts the situation in the 3D case, where a considerable amount of pioneering work in-volving quadrupolar particles at finite temperatures [15][16][17][18][19] and in the ground state 20 is available. On the other hand, there are a number of simulation studies dealing with non-quadrupolar, shape-anisotropic particles in 2D, examples being hard ellipsoids, 21-24 rectangles 25 or hard spherocylinders. 26 In particular, the density-driven nematic-isotropic transition of hard ellipsoids was already investigated in Ref. 24 and turned out being continuous or of first order depending on the aspect ratio. We also mention a density functional study 27 of 2D systems of anisotropic particles exhibiting isotropic phases and phases with simultaneous orientational and translational order. Thus, the important open question to be explored in the present article concerns the interplay of shape anisotropy and quadrupolar intermolecular interactions.
To this end we employ two types of numerical calculations, that is, optimization techniques based on ideas of evolutionary algorithms (EA) and constant-pressure MD simulations, revealing the system's behavior in the ground state and at finite temperatures. To facilitate the calculations, the non-electrostatic part of the interaction between two molecules is modeled via a purely repulsive, relatively stiff (but not infinitely hard) pair potential that foots on the ellipse-like contact distance suggested by Berne and Pechukas. 28 The analytical form of this pair potential resembles an anisotropic soft-sphere potential, but has a higher exponent. 29 This potential turned out to be suitable for both type of numerical calculations. Van der Waals interactions between the molecules are entirely neglected. Although this may seem somewhat unrealistic for real molecules, there are several advantages of the resulting (repulsive) potential: it reduces the number of parameters, it allows (as we will demonstrate) for a transferability between results at different system parameters, and it makes it possible to establish a connection of the present study to investigations of anisotropic colloidal particles in confined geometry. [30][31][32] The remainder of this article is organized as follows. In Sec. II we introduce the model and the pair interactions in detail. Section III is devoted to the methods used to calculate equilibrium structures (at T = 0 and T > 0) and corresponding data analysis methods. The resulting ground state structures at different pressures and aspect ratios are presented in Sec. IV A, and the corresponding finite-temperature results are discussed in Sec. IV B. Finally, we summarize our findings in Sec. V. Appendices A to C provide additional, complementary information.

II. MODEL
In our investigations we consider a two-dimensional system, embedded in the (x, y)-plane, where the anisotropic particles are only allowed to rotate around the z-axis. We use reduced units throughout, introducing parameters σ 0 for the length scale, 0 for the energy scale, and m 0 for the unit of mass.
Our particles are assumed to have an elliptic shape, characterized by the lengths of the two main axes, σ and σ ⊥ ; the indices refer to the orientation of the corresponding axes relative to the linear quadrupolar moment to be introduced below (see Figure 1). Defining a dimensionless shape anisotropy parameter, κ, we impose via that the surface area of our particles is independent of the actual value of κ. The isometric case of a disc is recovered for κ = 1.
The interaction between two-particles (with indices i and j) can be split into a short-ranged, anisotropic repulsion, V sr (r ij ,û i ,û j ), and a long-ranged potential, V lr (r ij ,û i ,û j ); the latter one stems from the interaction of the linear quadrupoles that the particles carry. We define the distance r ij = |r i − r j | between the centers of the particles and furtherr ij = r ij /r ij . Theû i and u j are the normalized orientation vectors of the respective linear quadrupolar moments. Within our model the shape of the particle is the same for κ and 1/κ, while the respective orientations of the embedded quadrupole moment are perpendicular for the cases κ and 1/κ (see bottom panels of Figure 1). Consequently, by choosing values of κ both smaller and larger than unity we are able to consider these two orientations of the quadrupolar moment within the same model.
With the above vectors at hand we introduce For the short-range contribution of the inter-particle interaction (subscript "sr") we use an anisotropic, repulsive potential: it has the simple functional form of an inverse power law (IPL) interaction, while its dependence on the connecting vector, r ij , and on the orientations of the quadrupolar moments,û i andû j , is inspired by a Gay-Berne potential: 34 The electrostatic part of the inter-particle interaction is based on linear quadrupole moments, for which the quadrupole tensor, 35Q , becomes in its eigenbasiŝ with Q being the strength of the moment. The interaction between two linear quadrupole moments can be written as 16,35 V lr (Q 2 , r ij ,û i ,û j ) = 1 4πε pm 3 4 with a, b, and c defined above and ε pm being the vacuum permittivity. Even though the interaction decays slower with the distance (i.e., ∼ 1/r 5 ) than, e.g., van der Waals interactions (∼ 1/r 6 ) we can still calculate inter-particle energies with sufficient accuracy via a real-space lattice sum. Here we have used a cutoff radius R cut /σ 0 = 30 for both, electrostatic and repulsive interactions. The interaction of two linear quadrupoles depends on their respective orientations. In contrast to dipoles, which tend to line up head-to-tail, isolated linear quadrupoles prefer the so-called T-configuration for κ close to unity (see Figure 2). A more detailed discussion of the impact of shape anisotropy on the two-particle arrangement with a minimal electrostatic energy is given in Appendix A.
Throughout, our numerical calculations were performed in the N P T -ensemble, where the relevant thermodynamic potential for vanishing temperature is the 1 FIG. 2. Typical arrangement of isolated quadrupolar particles with weakly anisometric shape (i.e., κ close to unity), induced by the fact that neighboring particles tend to arrange in a mutually orthogonal orientation. enthalpy, H, here E is the internal energy, P is the pressure, and S 0 is the area occupied by the system. Further, we introduce the temperature of the system, T . Based on the scales of length (σ 0 ), mass (m 0 ), and energy ( 0 ), the quadrupole strength, Q 2 , the pressure P , the temperature T and the time t are expressed via their respective reduced units, (Q * ) 2 = Q 2 /(4πε pm σ 5 0 0 ), P * = P σ 2 0 / 0 , T * = k B T / 0 (with k B being the Boltzmann constant) and t * = t/ σ 2 0 m 0 / 0 .

III. METHODS AND DATA ANALYSIS
In our approach we investigate in a first step the selfassembly scenarios of our particles at zero temperature (i.e., the ground state configurations). This information helps us to classify archetypical particle arrangements (each characterized by a particular spatial and orientational order), which, in turn, help us to understand the strategy of the particles of how to assemble in the energetically most favorable manner. For these investigations we employ an optimization tool based on evolutionary algorithms, detailed below.
In a second step we take advantage of these ground state configurations and use them as starting configurations in molecular dynamics (MD) simulations, performed at small, but finite temperatures. These investigations provide information on the thermodynamic stability of the ordered phases.

A. Evolutionary algorithms
In our effort to identify ground state configurations of our systems, we use an optimization tool based on ideas of evolutionary algorithms (EAs). 36 EAs are heuristic approaches designed to search for global minima in high-dimensional search spaces and for problems that are characterized by rugged energy landscapes.
In an effort to be compatible with the requirements of an N P T -ensemble, we introduce in our approach a unit cell of variable area and shape which creates (together with its periodic images) a system of infinite extent. In the desired configuration (the so-called ground state configuration) particles are located and oriented in this cell in such a way as to minimize the internal energy of the system, which at vanishing temperature is equivalent to the enthalpy.
We initialize the algorithm by creating a set of configurations where particles are located in the cell at random positions and have random orientations. These arrangements are graded by their respective fitness value, a quantity that provides evidence on how suitable this configuration is to solve the optimization problem. Since we are interested in finding ground state structures, a high fitness value of a particular configuration corresponds to a low value of the enthalpy per particle.
Then we iteratively use existing particle arrangements to create new ones by applying one of the two following operations: crossover and mutation. In the former one we first select two configurations where the choice is biased by high fitness values of the two configurations. Characteristic features of both particle arrangements (such as lattice vectors, particle positions, and particle orientations) are then combined to form a new configuration. The mutation operation, on the other hand, introduces random changes to a randomly chosen configuration, such as moving or rotating an arbitrarily chosen particle or changing the lattice vectors. Typically 2000 iterations are required for a particular state point until convergence towards the global minimum has been achieved.
Our implementation of EAs is memetic, i.e., we combine global and local search techniques: each time a new configuration has been created with one of the two above mentioned EA operations, we apply the L-BFGS-B 37 algorithm which guides us to the nearest local minimum. This algorithm has been applied successfully for a broad spectrum of systems, both in two and three dimensions (see, for instance, Refs. [38][39][40][41][42][43]. Within the context of the present contribution, it should also be mentioned, that a suitable extension of the formalism is also able to cope with long-ranged interactions, as they are encountered in charged systems (see Refs. 44 and 45).
For the present contribution we have performed with our EA approach computations for 211 evenly-spaced values of κ ∈ [0.4, 2.5] and for several combinations of (Q * ) 2 ∈ {0.2, 2, 20} and P * ∈ {0.1, 1, 10}. Due to computational limitations we have considered unit cells that contain up to twelve particles.

B. Molecular dynamics simulations
In order to investigate structural changes of our system at finite temperature, we perform molecular dynamics (MD) simulations at constant pressure and temperature for an ensemble of N = 840 particles. To this end we use the Berendsen weak-coupling scheme, 46 where the classical Newtonian equations of motion are supplemented by terms representing the coupling to a heat-and a pressurebath.
In the following, we first introduce the translational equations of motion proposed by Berendsen et al., 46 specializing to the case of a two-dimensional system. Considering a particle with index i (i = 1, . . . , N ) one haṡ where the vectors b α (α = 1, 2) define the simulation box and its shape.
In the above equations, the dot represents a timederivative. Further, r i and v i denote the position and the velocity of particle with index i, respectively; F i is the force on particle i with mass m 0 , and K is the compressibility. The pressure tensors, P (actual pressure) and P (target pressure), are defined below. Further, T trans is the actual kinetic temperature, where the denominator 2(N −1) represents the 2N translational degrees of freedom minus two constraints due to momentum conservation in each spatial dimension. To ensure this conservation during the numerical integration, we set the total momentum to zero every 100 time steps.
The (time) constants τ trans and τ P determine the speed of relaxation of T trans and P towards their respective target values defined by the bath. Specifically, the corresponding coupling equations reaḋ where the pressure tensors are defined as and Here, ⊗ denotes the dyadic product of two vectors, S 0 stands for the area of the simulation cell, 1 is the unittensor, and F ij is the force on particle i exerted by particle j.
To describe the rotational motion of the particles, Eqs. (9)- (14) are supplemented with the following differential equations for the angular velocity, ω i , of particle i:ü Here, M i and I are the torque on particle i and the corresponding moment of inertia of particle i, respectively; × denotes the vector product. The kinetic temperature for the rotation is defined as k B T rot = i Iω 2 i /N . This temperature is controlled in analogy to Eqs. (12).
To solve the translational equations of motion, i.e. Eqs. (9)-(14) we use a modified leap-frog integrator as proposed in Ref. 46. For the solution of the rotational equations of motion, Eqs. (15)- (16), an analogous procedure is performed (see Ref. 47).
To initialize the particle positions in our system, we use the unit cells obtained from the ground state calculations in the preceding investigations (see Subsec. III A). Specifically, we take for each (κ-dependent) ground state the respective unit cell and arrange copies of this cell such that the simulation cell has a minimal circumference. For this simulation cell (which is now defined by the simulation box vectors b 1 and b 2 ) periodic boundary conditions are applied. The pair forces F ij and torques M ij are truncated at R cut = 6σ 0 and are then shifted in order to make them vanish smoothly.
The MD calculations have been performed at several values of the reduced temperature T * , that is, T * = 0.1, 0.2, . . . , 1.6. In our simulations we used a time increment ∆t * = 0.00136172 (corresponding at a macroscopic level to ∆t ≈ 2 fs). An adequate choice of the time step is imposed by the mass of the particles, for which we have assumed the mass of benzene, m 0 = 78u. Each simulation extends over 210 000 time steps. During the first 20 000 MD steps, the target temperature is gradually increased from T * = 0 towards the respective target value. Structural quantities are then extracted only during the final 10 000 time steps of the simulation. The values for the temperature-and pressure-coupling constants τ trans , τ rot and τ P were chosen to be 80 time steps. Of course, the actual volume change of the box is also influenced by the compressibility; its reduced, dimensionless counterpart, K * = Km 2 0 σ 2 0 / 0 , is set to 0.01. We assume the mass distribution within the particles to be homogeneous and thus obtain for the dimensionless moment of inertia I * = 0.25 κ + 1 κ .

C. Structural analysis -order parameters
In order to quantify both the positional and the orientational order of the particles we introduce three different sets of order parameters: (i) The positional order can be quantified via twodimensional bond-orientational order parameters (BOOPs), 48 introducing, in addition, weight factors 49 that are related to the side lengths of the Voronoi polygon around each particle: (17) here (j ∈ N i ) indicates that particle j is a nearest neighbor of particle i and l ij is the length of the side of the Voronoi polygon that separates the two particles. φ ij is the angle enclosed by the bond between particles i and j and the reference axis, which we choose as the x-axis. For the parameter n we have chosen the values 4 and 6, highlighting thus via Ψ 4 or Ψ 6 four-or six-fold symmetry, respectively. While BOOPs are calculated for only the best configuration at vanishing temperature, averages over several configurations are taken in MD simulations, indicated by the brackets in Eq. (17) and in the following.
(ii) An obvious candidate for quantifying the orientational order of the particles is the nematic order parameter S, which is introduced via the tensor order parameter T, defined as S is then simply the positive eigenvalue of the tracless tensor T, which can easily be calculated as The eigenvector associated with S is called the directord and indicates a preferred orientation in the configuration. Note that S is always positive, except in the case that all elements of T vanish: then no preferred direction is present in the system.
A different way of how to measure the orientational order is via the parameter β, defined as (iii) In addition, a variety of order parameters combining positional and orientational order can be defined. An example of such an order parameter, which has been used in a previous contribution 50 is Concluding, we emphasize that orientational order refers to the orientation of the particles; the orientation of the linear quadrupolar moment is then imposed by the value of κ: if κ > 1, then the orientation of the moment is parallel to the main axis of the particles, while for κ < 1 these two orientations are perpendicular to each other.

A. Ground state configurations
In our discussion of the ground state configurations we first present the different structural archetypes that we could identify with our EA-approach. We then present the diagram of states: it provides information on the κ-ranges where the respective structures are the energetically most stable ones; we further discuss how the thermodynamic properties and the order parameters of our systems vary over a representative range of κ.

Structural archetypes
a. Structures with one preferred orientation For κvalues both considerably larger or smaller than unity (i.e., by a factor of ≈ 2), the particles tend to align parallel to each other in case of strong anisotropy: they accomplish this by forming tilted rows, characterized by a high density of particles along these lines; these arrangements are denoted as parallel displaced-configurations ("PD"; see Figure 3). Since quadrupolar particles avoid a headto-tail arrangement, neighboring rows repel each other, inducing thereby large gaps between these lanes; the parallel offset between neighboring rows can be very sensitive to small changes in κ, leading to a slight modulation of the BOOPs as functions of κ (to be discussed below). An interesting special case of this structure is observed for a vanishing quadrupolar moment (i.e., (Q * ) 2 = 0), where particles form a distorted hexagonal lattice, denoted as the parallel-configuration ("P"; see top left panel of Figure 3).
b. Structures with two preferred particle orientations Here the preferred structural arrangements are throughout herringbone-configurations ("HB") where particles form alternating, parallel lanes, each of them being characterized by a specific particle orientation: in Figure  4 particles pertaining to different lanes are colored red and blue, respectively. The relative orientation between neighboring particles depends in a highly sensitive manner on κ.
Apart from two special cases (that are encountered for κ-values very close to unity and which will be discussed below) two different versions of the HB-configurations emerge from a closer analysis of the obtained structures: (i) for κ-values somewhat closer to unity, we observe a rather dense structure, which we denote HB dense -  configuration (see panels in the central column of Figure 4); (ii) for κ-values that differ more strongly from unity, particles arrange in densely-populated rows of alternating orientation (see panels in the right column of Figure 4). Since neighboring rows repel each other, this structure has a lower overall den-sity which we therefore denote as HB loose .
The two special cases of the HB structure mentioned above are observed for κ-values very close to unity where neighboring particles prefer strict mutual orthogonal orientations with respect to their nearest neighbors. (i) For large (Q * ) 2 -and small P * -values, a perfect arrangement of mutually orthogonally oriented particles with an un-derlying square pattern can be observed, denoted as the square T-configuration ("T sq "; see left panel of Figure  4). (ii) Further, a closely related arrangement has been identified which is now based on an underlying hexagonal pattern; it is denoted as the hexagonal T-configuration ("T hex "; see bottom left panel of Figure 4); for this particular configuration a mutually perfect orthogonal orientation of nearest neighbors can only be realized for κ = 1.
c. Structures with three preferred particle orientations A very interesting phenomenon observed in our system is the occurrence of structures where particles arrange in three preferred orientations, characterized by a vanishing nematic order parameter S.
The spatial particle arrangement is here -irrespective of the values of Q * and P * -reminiscent of a trihexagonal tiling, 51 consisting of regular hexagons that are connected by triangles; we thus denote it as the "TH"-configuration (see Figure 5 where also the hexagons and the triangles are highlighted). Throughout, the relative angle between the orientations of neighboring particles is π/3, reflected in the order parameter which assumes a value of β = cos 2 (π/3) = 1/4 (see Figure 7).
We note that the occurrence of the TH-configuration strongly depends on the choice of (Q * ) 2 and P * ; in some cases, this particle arrangement is not observed at all (see Figure 8). In contrast, variants of the HB-configuration can almost always be observed in some range of κ (see discussion below).
d. More complicated structures Using the EAapproach we also identify more complicated structures, not conforming to the mechanisms described above. These structures require a larger number of particles per cell and turn out to be stable only within very small ranges of κ (see Subsec. IV A 2). While we observe several different variants of such particle configurations, they share a common feature, namely a mesh-like lattice, where chains of particles undulate back and forth (see Figure 6). We denote them as branched structures ("B").

Diagram of states
We now discuss the diagram of states which summarizes the occurrence of the previously identified archetypes of ground state configurations of our system for selected values of (Q * ) 2 and P * as we vary κ. Since most of the interesting features of our systems can already be captured in the diagram of states obtained for (Q * ) 2 = 2 and P * = 1, we focus from now onwards on this particular set of parameters: we present data for the enthalpy and for a selection of appropriate order parameters introduced in Subsec. III C, that help to identify the respective structures. The filling fraction η, which is also displayed in the following figures, is defined as η = N πσ 2 0 /4/S 0 . At this point we remind the reader that the shape (and thus the orientation of the particles) is invariant under the transformation κ ↔ 1/κ while the relative orientations of the quadrupolar moments for the cases κ and 1/κ are mutually orthogonal: with this symmetry in mind, we can easily disentangle the impact of (Q * ) 2 and P * on the structure formation by comparing the relevant physical properties for the cases κ and 1/κ. These properties are displayed for (Q * ) 2 = 2 and P * = 1 in Figure  7 along with a color-coded, horizontal bar (above this panel) which indicates these κ-ranges where the respective particle configurations are the energetically most favorable ones.
On a qualitative level we observe that the enthalpy curve is continuous over the entire investigated κ-range; however, it shows kinks at particular values of the aspect ratio which provide a first evidence for the occurrence of discontinuous transitions between the ground state configurations. The specific enthalpy, H * /N , shows a pronounced local minimum at κ = 1 (i.e., for circular particle shapes): here the T-configuration is dominant, verified by the fact that in this κ-region Ψ 4 = 1 and Ψ 6 = 0. In addition, we observe S = 0 and β = 0 = cos 2 (π/2), indicating thereby a relative orthogonal orientation of neighboring particles. In contrast, for more anisometric particle shapes (i.e., small and large κ-values), the enthalpy rapidly decays to rather small values; actually H tends to minus infinity for κ → 0 and κ → ∞ as the quadrupoles start to overlap and the repulsive soft core shrinks as σ ⊥ → 0. For small and large κ-values the formation of parallel rows is energetically most favorable and PD-configurations are observed. Here the nematic order parameter S is the appropriate quantity to characterize the emerging structure: indeed, S assumes in these κ-regions the value 1, indicating that a single orientation prevails.
In contrast, for the two intermediate κ-ranges, where the enthalpy assumes local maxima the ground state configurations for κ and 1/κ are distinctively different; for these κ-values we observe the formation of more complex structures, reflected by a rather intricate variation of the different order parameters with κ: (i) Increasing first κ beyond the region where the T-configuration is stable, the system changes -after a very narrow κ-region where the HB dense -structure occurs -via a discontinuous structural transition (identified by a jump in η) into the TH-configuration; the latter one is characterized by S = 0, Ψ 6 = 1 and β = 0.25 = cos 2 (π/3) which provides evidence that the difference in orientation angles between nearest neighbors is π/3. Upon further increasing κ we pass again a very narrow interval where a branched structure is stable and we eventually reach -again via a discontinuous structural transition -the HB loose -structure; for this configuration none of the order parameters assume any characteristic value. Eventually we identify for even larger κ-values the aforementioned PD-structure.
(ii) Decreasing, on the other hand, the value of κ beyond the range where the T-structure is stable, we observe HBconfigurations: first the one with the larger density (i.e.,  the HB dense -structure) and then the HB loose -structure. The shape of β (see Figure 7) indicates that the transition between HB dense and HB loose (and also the subsequent transition to the PD-structure) is of first order, i.e., at some point the relative orientation of neighboring particles is discontinuous.
At this point it should be mentioned that η and the order parameters do not necessarily behave in the same manner as the system changes from one structure to the other: some order parameters or η may change abruptly, indicating a first order phase transition, while the other parameters change continuously. As an example we refer the reader to the transition HB loose → PD for κ ∼ 0.45.
Another interesting feature emerges as we compare the curves of the order parameters shown in Figure 7 obtained for the parameters (Q * ) 2 , P * = (2, 1), with the corresponding data calculated for the parameter sets (Q * ) 2 , P * = (0.2, 0.1) and (Q * ) 2 , P * = (20, 10), shown in Figures 8 and 9, respectively. The corresponding curves reveal striking similarities, suggesting that appropriate scaling relations of the order parameters via the values of the quadrupole moment and the pressure hold for the respective ground states. In contrast (and interestingly), the enthalpy curves obtained for the different sets of data differ rather substantially in magnitude. We will discuss a possible background scenario of these observations in more detail in Appendix B. There we will show that indeed a scaling relation for the enthalpy of a systems of hard, quadrupolar particles can be derived (see Appendix B 1). However, it seems that the softness of our particles -remember that we consider in this contribution particles with a soft (albeit rather steep) coreleads to a breakdown of this scaling law. This feature is presumably due to the fact that a simultaneous scaling of the quadrupole moment and of the pressure also induces a change of the density. To take into account this effect properly, we suggest in Appendix C 1 an empirical scaling law for the ground state enthalpy of a system of soft ellipsoidal particles and provide numerical evidence for its justification.

B. Results at finite temperatures
In the current subsection, we analyze the MD results that we have obtained for all considered aspect ratios κ and a set of reduced temperatures T * . As previously mentioned, the ground state configurations that we have obtained with the help of the EA algorithm for a variety of aspect ratios, serve now as initial configurations for the MD simulations. Corresponding snapshots for the most interesting and most representative configurations, obtained after equilibrating the system, are presented in Figure 10.
Before starting the investigations, it is worth to briefly consider typical values of (Q * ) 2 , P * in experimental systems. As an example, we consider the quadrupolar Gay-Berne model of benzene proposed by Golubkov et al. 13 . We use their quadrupole strength, Q benzene = −30.5812 · 10 −40 Cm 2 , and their spherical diameter, σ 0 ≈ 0.307 nm (the latter value is based on the Gay-Berne contact distance). Using the surface tension of water (the 2D equivalent of pressure) defined in Ref. 52, P water = 71.99 10 −3 Nm −1 at 25 • C, we arrive at the ra-tio (Q * benzene ) 2 /P * water = Q 2 benzene /(4πε pm σ 4 0 P water ) ≈ 4.5.
However, since the point quadrupole approximation is known to overestimate the interaction strength for narrow interparticle configurations 13,53,54 due to the singularity in V lr [see Eq. (7)], we consider here a reduced value of (Q * ) 2 /P * = 2 instead of 4.5. With this choice we hope to cover not only benzene molecules, but also other quadrupolar molecules. By identifying P * = 1 with P water , we arrive at the following energy scale 0 = P water σ 2 0 /P * = P water σ 2 0 = 4.09 kJ/mol. In the following, we present MD simulation results for the parameter pairs (Q * ) 2 , P * = (2, 1), and (4, 2). We start with a detailed investigation of the structure at (Q * ) 2 = 2, P * = 1 and different temperatures.
The MD simulations reveal that the structures predicted for vanishing temperature remain stable also at low temperature (T * = k B T / 0 = 0.2, see panels in the second column of Figure 10). As T * increases monotonously, defects start to form (see panels in the third column of Figure 10 a slight wave-like modulation of previously straight lines. Finally, once the temperature has been raised above a certain threshold value (see panels in the fourth column of Figure 10), the crystalline order is rapidly lost. The corresponding transition temperatures depend strongly on κ: configurations with κ ≈ 1 (T sq -configuration) and κ far from unity (PD-configuration) turn out to be the most stable ones: in this situation, particles can approach very closely and exert thus a strongly attractive quadrupole-quadrupole interaction; these ordered structures break up only at temperatures as high as T * = 1.1 and T * = 1.2, respectively. In contrast, the rather complicated B-configurations melt already at temperatures as low as T * = 0.5 (with κ = 1.66).
In order to investigate and to locate the transition of the system from the ordered to the disordered regime, we focus in the following on the reduced potential energy, E * pot , and the reduced system area, S * 0 , (in units of 0 and σ 0 , respectively). Figure 11 depicts E * pot and S * 0 as functions of the reduced temperature T * for all considered values of the aspect ratio κ. Upon increasing the temperature, E * pot progressively decreases in magnitude and finally approaches zero, reflecting the diminishing role of particle interactions. At the same time, the area S * 0 increases. Interestingly, we observe along this process that both E * pot (T * ) and S * 0 (T * ) show discontinuous changes within very small temperature intervals for all κ-values investigated; these "jumps" are found for both quantities at approximately the same temperature. We attribute these observations to the occurrence of a first order phase transition. The temperatures that delimit these intervals are marked in Figure 11 by contour lines as functions of κ. Only for κ 1.66 no such discontinuity in E * pot (T * ) could be resolved; this fact can be related to the finite size of the temperature grid. Intentionally, we have not evaluated the susceptibility, since it is well established that energy fluctuations are not correctly reproduced within the Berendsen scheme. 55, 56 We now arrive at the discussion of the previously defined BOOPs Ψ 4 and Ψ 6 [see Eq. (17)]; they are displayed in Figure 12. As already shown in Figure 7, we find for all our ground state configurations intervals in κ where at least one of the two parameters does not vanish. Considering now the temperature dependence of these quantities, we observe for all κ-values investigated a discontinuous change of both Ψ 4 and Ψ 6 (if not zero in the ground state) at exactly the same temperatures where a discontinuous change in S * 0 was observed. To analyze the structure of the system at finite temperatures on an even more quantitative level, we have calculated in addition various coefficients of the pair distribution function in an expansion in terms of rotational invariants. To be more specific, we use a two-dimensional version of the threedimensional coefficient functions of the pair distribution function, 57 i.e.
Here we focus on the coefficients g 000 (R), g 220 (R), and g 202 (R) defined via the following two-dimensional rotational invariants: The function g 000 (R) corresponds to the familiar pair correlation function. The other functions, g 202 (R) and g 220 (R), describe, in addition, a local orientational order of the particles. Similar to 3D systems 58 we can interpret these coefficients as follows: g 220 (R) provides information about the conditional probability density of a particle (relative to the bulk probability density) whose orientation axis is aligned in parallel (positive value) or orthogonal (negative value) to the orientation axis of a con- sidered particle; g 202 (R) describes the conditional probability density of a particle (again, relative to the bulk probability density) positioned along or aside the axis of a considered particle.
In the panels of Figure 13, we present numerical results for these three functions for κ = 2.1 (see also the corresponding snapshots shown in the panels in the bottom row of Figure 10). For all three correlation functions the first peak is located at around 0.8 σ 0 : this position does not exactly mark the face-to-face alignment of the particles (which occurs at ∼ 0.69 σ 0 ) and thus provides evidence for a slightly shifted parallel configuration, induced by the quadrupole (see snapshots shown in the panels of the bottom row in Figure 10). For T * = 1.2, we observe a rapid decay of g 000 (R) towards unity for large particle distances R: the crystalline order has completely van- ished, reflected by the missing peaks for larger R-values.
Since at low temperatures all particles are oriented in the same direction, g 220 (R) is equivalent to g 000 (R); however, as the temperature attains T * = 1.2, we observe that for large R-values the orientational order is lost and thus g 220 (R) vanishes already at R 2σ 0 . Finally, g 202 (R) oscillates around zero in the ordered phase for low tem-peratures (T * < 1.2) due to the fact that g 202 (R) isby definition -not able to contribute to the overall particle density. For T * 1.2, g 202 (R) completely vanishes for large R-values (R 2σ 0 ): obviously, a loss of orientational order occurs for T * 1.2. This temperature threshold, which can be interpreted as a melting temperature, agrees fairly well with the temperature where the discontinuous changes in E * pot (T * ) and S * 0 (T * ) (see Figure 11) and in the BOOPs (see Figure 12) were observed. Concluding, we note that related investigations carried out for other κ-values led to analogous conclusions about the corresponding melting temperature. These observations motivate investigations on a melting curve that separates the ordered from the disordered phase as we increase the temperature. This line is displayed in Figure 14 for all considered aspect ratios κ. We emphasize that these data represent only an estimate for the true, two-phase coexistence lines characterizing a first-order transition. For the latter, one would also expect the occurrence of a hysteresis, i.e., the observation of two different curves, depending on whether the system is heated up or cooled down from a low or a hightemperature state, respectively. We briefly come back to this issue below. Another interesting feature of the melting curve that can be observed is that it exhibits some degree of symmetry in shape when exchanging κ and 1/κ (see Figure 14). Of course, we would not expect full symmetry since the electrostatic properties for the cases κ and 1/κ are different (see bottom panels in Estimate of the melting curve T * melting ((Q * ) 2 =2, P * = 1, κ) that separates the ordered phase (at low temperatures) from the disordered phase (at high temperatures); for details cf. text. Data along on the line correspond to results obtained for disordered state points. Figure 1). From the data we can conclude that for particles with large eccentricities the melting occurs at higher temperatures than for κ-values close to unity.
Given the large amount of numerical calculations required to construct the melting line in Figure 14, which was calculated for one particular set of parameters, (P * , (Q * ) 2 ), it would be obviously desirable and helpful to have a scaling relation at hand which allows to easily obtain (or to extrapolate) corresponding melting lines for other parameter sets. In Appendix B 2 we show that such a relation does indeed exist for hard particles. This relation states that the probability to encounter a microscopic configuration of the many-particle system in phase space is invariant under the simultaneous transformations Q 2 → µQ 2 , P → µP , T → µT , with µ being a scaling factor. However, in this contribution we consider soft particles (even though characterized by a rather harsh repulsion, see Eq. (2)). Nevertheless, as we discuss in Appendix C 2, it is also possible for the system at hand to provide via a suitably adapted scaling law a rough estimate of the location of the melting line at scaled parameters.
Finally and for the sake of completeness, we now discuss our investigations on the melting transition as obtained in a cooling process ("simulated annealing"). The central question that we address is whether the ground state structures can be reproduced -at least locallywith MD simulations starting from a disordered phase at higher temperatures. To this end we performed simulations at (Q * ) 2 = 2 and several values of κ, initializing the system with random particle positions and orientations and cooling it down gradually. The initial pressure and temperature were set to P * = 10 and T * = 5, respectively; the initial box-shape is quadratic with a side length of 50σ 0 . Within the first 200 000 MD steps we linearly decreased the pressure and the temperature down to P * = 1 and T * = 0.1, respectively. The simulations extended in total over 410 000 MD steps. For the other simulation parameters we refer the reader to Sec. III B. Our data provide evidence that the predicted ground state could be obtained via this simulated annealing process only for very few state points. In general, the formation of ordered structures via such a process is hampered and delayed by frustration effects, especially for κ-values far from unity, where particles encounter -due to their elongated shapes -difficulties to rotate. As an example we discuss the pair correlation function g 000 (R), obtained for κ = 0.8 (i.e., a value close to unity) and depicted in Figure 15(a). We observe a coincidence in the peak positions of g 000 (R), but not in their heights. We interpret this deviation as an artefact of the crystallite structure appearing after cooling (see Figure 15(b) and (c)). This might be a consequence of the lack of long-range order in 2D systems.

V. CONCLUSIONS
In summary, we have proposed a simple model that mimics the essential features of elongated, organic molecules without a net charge or dipole moment: it consists of soft ellipsoids with an embedded linear quadrupole moment. The preferred orthogonal arrangement of linear quadrupole moments in close proximity to each other represents an interesting contrast as compared to dipolar systems, where particles prefer paral-lel arrangements, thus often forming chains. The selfassembly scenarios of our system result from a competition between the shape anisotropy of the ellipsoids and the quadrupolar interactions.
Operating in the NPT ensemble, we have investigated the system for numerous different sets of system parameters, (Q * ) 2 and P * , and could identify different strategies of the system, depending on the shape anisotropy κ. While we always observed configurations of parallel rows of particles for κ-values far from unity, the sequence of structures between these two limiting cases strongly depends on the competition between the short-ranged and the long-ranged, electrostatic contributions to the interaction. For κ ≈ 1 we observe two different configurations with orthogonal particle arrangement, one of them based on a square lattice, the other closely related to a hexagonal lattice. Intermediate ranges of κ -both for κ > 1 as well as for κ < 1 -are mostly dominated by two variations of the herringbone structure. In addition, we observe for selected κ-values a non-trivial lattice, closely related to the trihexagonal tiling. In very small ranges of κ, more complicated, branched structures can emerge. However, these turn out to be in general rather unstable at finite temperatures, as shown in complementary MD simulations: thus we speculate that they are metastable at vanishing temperature.
For future work, it would be useful to extend our simple model in two directions: (i) the particles interaction should include a general quadrupole moment and (ii) the system could be confined to a slab geometry or be extended to full three dimensions. These extensions would nicely meet the considerable experimental interest in the self-assembly of complex organic molecules on surfaces. 59,60 Some of these systems show a remarkable variety of different structures, which can even be controlled via external parameters such as magnetic fields 61 or light with different polarization, 62 opening thereby the route to many interesting technological applications. [63][64][65][66][67][68] The understanding of our simple model considered in the present contribution could serve as a suitable starting point for related investigations of the more complex molecules considered in these studies. Extending the model step-by-step towards more complicated shapes and to more intricate effective interactions would help to understand the numerous competing effects. 5,69 Finally we note that so-called Inverse Patchy Colloids (IPCs), 70 i.e., colloids decorated with charged patches, represent due to their charge distribution a closely related system, as they also carry a linear quadrupole moment. Appendix A: Two particle ground state configurations at particle contact We consider two quadrupolar particles at contact, i.e., r 12 = σ(r 12 ,û 1 ,û 2 ) [see Eq. (4)] and optimize their respective and relative orientations by minimizing their electrostatic energies. In Figure 16 we display the angles α 1 and α 2 , enclosed by the vector connecting the centers of the two particles and the orientations of the respective quadrupolar moments in the corresponding optimal configurations, as functions of κ (for the definitions of the angles see also inset in Figure 16). Obviously, the Tconfiguration is only stable for values of κ close to unity. For strong anisometry (κ 0.75 or κ 1.2) the ground state configuration is parallel displaced (PD).

Appendix B: Reduction of parameter space for hard particles
The goal of the present subsection is to derive scaling relations for the system parameters Q 2 (quadrupole strength) and P (pressure) under which the thermodynamic properties of the system (and thus its ground state configurations) remain unchanged. Of course, these relations also hold for our reduced parameters, which we introduced at the end of Sec. II.
We assume in the following a simplified version of our interaction potential, namely the quadrupolar hardellipse (QHE) model, defined by its interaction σ(r ij ,û i ,û j ) has been defined in Eq. (4). Thus, contrary to our original system (see Sec. II), the QHE model consists of impenetrable ellipsoidal particles, whereas the quadrupolar, long-range part of the interaction is identical to our original model, Eq. (7). The potential energy, E pot , of a specific configuration of the QHE system is defined via where r N = (r 1 , ..., r N ) andû N = (û 1 , ...,û N ) specify the microscopic configuration of the particles.

Ground state (vanishing temperature)
At vanishing temperature, the particles are positioned and oriented in a non-degenerate ground state in an NPT ensemble such that the enthalpy, H, is minimal, that is The scaling properties of the QHE interaction potential (via Eq. (B1)), with the quadrupole strength Q 2 induces a scaling property of the entire potential energy [cf. Eq. (B2)], i.e., E pot (µQ 2 ) = µE pot (Q 2 ); henceforward, µ is a simple scaling parameter. Since the enthalpy H is a linear combination of E pot and P , replacing Q 2 by µQ 2 [cf. Eq. (B3)] leads to the same ground state only if the pressure also scales with µ. Thus, we can conclude that a particular ground state configuration of the QHE model (specified by {r N GS ,û N GS , S GS 0 }), remains invariant under the transformations Q 2 → µQ 2 and P → µP .

Finite temperatures
We next consider the QHE system at T > 0. The phase space probability for the occurrence of a microscopic configuration (specified by {r N , u N }) in a phase space volume of infinitely small extent in an NPT-ensemble is given by where Z conf is the corresponding partition function of the ensemble. It is obvious that Eq. (B4) is invariant under the transformations E pot → µE pot , P → µP , and T → µT . Since we know from Eqs. (7) and (B2) that the potential energy of non-overlapping, impenetrable particles is proportional to Q 2 , i.e., we obtain the simple rule, that ρ(r N , u N , S 0 ) is invariant under the scaling law Q 2 → µQ 2 , P → µP , and T → µT . Thus, if we have a reference system (specified by Q 2 , P , and T ) at hand, we can simply calculate the properties of another system (index "new") via the following procedure. Given the size of the corresponding molecule, σ new , its aspect ratio κ and the quadrupole strength, Q 2 new , it is thus possible to estimate the order-disorder transition temperature for the new system: we first calculate the reduced quadrupole strength via (Q * new ) 2 = Q 2 new /(4πε pm σ 5 new 0 ) from which we obtain the scaling factor µ from µ = (Q * new /Q * ) 2 . The order-disorder temperature and pressure of the new system then follow from the corresponding quantities of the reference system, T * melting and P , via T * melting,new = µT * melting and P * new = µP * .
Appendix C: Applicability of parameter scaling for soft particles We now explore to which extent we can apply the parameter scaling relations derived for the QHE system to the actual system of soft particles investigated in the main part of this paper (see Sec. II).

Ground state (vanishing temperature)
The total potential energy for a specific configuration of soft quadrupolar ellipsoids is defined through E pot (r N ,û N ; Q 2 ) = i,j;i<j V sr (r ij ,û i ,û j ) + V lr (r ij ,û i ,û j ; Q 2 ). (C1) As pointed out in the discussion of the ground states in Sec. IV A, the order parameters reveal as functions of κ strong similarities when one compares their values as obtained for the parameter sets (Q * ) 2 , P * = (0.2, 0.1) and (Q * ) 2 , P * = (20,10). This suggests that the soft particle system fulfills at T = 0 a similar scaling relation as the QHE system. However, the corresponding enthalpy curves match only in terms of shape, but not in terms of magnitude. We thus introduce a phenomenological correction to the simple scaling law of the enthalpy as developed in Appendix B 1. Specifically, we assume that a scaling of P and Q 2 with a factor µ results in a small rescaling of the ground state coordinates, i.e. r i → γr i , where γ attains values close to unity. We further assume that, when passing from one ground state, obtained for parameters (Q 2 , P ) and a filling fraction η, to a new ground system (index "new"), specified by parameters (µQ 2 , µP ) and a filling fraction η new , the rescaling factor γ is related to the change in density according to As a consequence, the ground state cell volume is scaled with γ 2 . Regarding the total potential energy, we assume that its scaling with respect to Q 2 and the rescaling of the positions follows the same law in the QHE system [see Eqs. (B2), (B1) and (7)], yielding E GS pot,new (γr ij ,û i ,û j ; µQ 2 ) ≈ µ γ 5 E GS pot (r ij ,û i ,û j ; Q 2 ).
Collecting all contributions, we obtain the following approximate expression for the enthalpy: We check the applicability of this phenomenological scaling rule by calculating, as an example, enthalpies of a new system from the values of the enthalpies obtained at (Q * ) 2 = 0.2, P * = 0.1 and (Q * ) 2 = 20, P * = 10.
The corresponding values for the enthalpy as functions of κ are presented in Figure 17. From these data we can conclude that Eq. (C4) seems applicable for κ-values close to unity where the corresponding structures possess a 4-fold symmetry (T-configuration). For particles with larger aspect ratios, substantial deviations between the respective enthalpy curves occur.

Melting curve
As shown in Appendix B 2, the properties of the QHE system fulfill a parameter scaling relation also at finite temperature.