Understanding the Formation of PbSe Honeycomb Superstructures by Dynamics Simulations

Using a coarse-grained molecular dynamics model, we simulate the self-assembly of PbSe nanocrystals (NCs) adsorbed at a flat fluid-fluid interface. The model includes all key forces involved: NC-NC short-range facet-specific attractive and repulsive interactions, entropic effects, and forces due to the NC adsorption at fluid-fluid interfaces. Realistic values are used for the input parameters regulating the various forces. The interface-adsorption parameters are estimated using a recently introduced sharp-interface numerical method which includes capillary deformation effects. We find that the final structure in which the NCs self-assemble is drastically affected by the input values of the parameters of our coarse-grained model. In particular, by slightly tuning just a few parameters of the model, we can induce NC self-assembly into either silicene-honeycomb superstructures, where all NCs have a f 111 g facet parallel to the fluid-fluid interface plane, or square superstructures, where all NCs have a f 100 g facet parallel to the interface plane. Both of these nanostructures have been observed experimentally. However, it is still not clear their formation mechanism, and, in particular, which are the factors directing the NC self-assembly into one or another structure. In this work, we identify and quantify such factors, showing illustrative assembled-phase diagrams obtained from our simulations. In addition, with our model, we can study the self-assembly dynamics, simulating how the NCs ’ structures evolve from few-NCs aggregates to gradually larger domains. For example, we observe linear chains, where all NCs have a f 110 g facet parallel to the interface plane as typical precursors of the square superstructure, and zigzag aggregates, where all NCs have a f 111 g facet parallel to the interface plane as typical precursors of the silicene-honeycomb superstructure. Both of these aggregates have also been observed experimentally. Finally, we show indications that our method can be applied to study defects of the obtained superstructures.

First experimentally observed a few years ago [73,85,86,88], the PbSe silicene-honeycomb 2D superlattice is expected to exhibit a Dirac-type band structure, with the semiconductor band gap preserved [89][90][91]. Hence, such a material would combine the properties of semiconductors with those of graphene, making it very interesting for optoelectronic applications. This superlattice is formed by PbSe NCs, approximately having a 5-nm size and a rhombicuboctahedron shape, disposed as in the silicene lattice, i.e., in a periodically buckled honeycomb lattice, and all oriented with a f111g facet parallel to the superlattice plane. The PbSe NCs, before the self-assembly, i.e., when still dispersed in solvent, are fully covered by oleic acid ligands, which are strongly attached at the f111g, f110g facets and weakly bonded at the f100g facets [92]. The formation process of the PbSe silicene-honeycomb superlattice typically occurs at the solvent-air interface. During the process, ligands desorb from the NC f100g facets, allowing oriented attachment between PbSe NCs by opposite f100g facets [73]. In similar experiments [73,85,87], self-assembly of NCs into square superlattices was observed, with all NCs oriented with a f100g facet parallel to the superlattice plane, and with oriented attachment between PbSe NCs by opposite f100g facets also performed. Although the relevant forces involved in the formation of these superlattices are known [93], the precise mechanisms are far from understood. In particular, it is not clear yet why the PbSe NCs form square superlattices in some cases and silicene-honeycomb superlattices in others. In this work, we shed light on this, introducing a new coarse-grained molecular dynamics (MD) model to study the self-assembly of NCs at fluid-fluid interfaces. By presenting simulations with our MD model, we illustrate how the intricate interplay between interface-adsorption forces orienting and keeping the NCs at the interface and attractive forces only between f100g facets of close-by NCs drive the self-assembly of PbSe NCs into either square or silicene-honeycomb superstructures. As a confirmation that our model correctly captures the key features of the NC self-assembly, in our simulations we observe, in the formation of both the silicene-honeycomb and square superstructures, the same few-NCs precursors observed experimentally, that is zigzag aggregates, where all NCs have a f111g facet parallel to the interface plane, and linear chains, where all NCs have a f110g facet parallel to the interface plane, respectively. In addition, the NC superstructures obtained with our MD model present the same kind of defects observed in the experimental superlattices, indicating that our MD model can also be used to study the formation and stability of defects in these structures.

II. COARSE-GRAINED MD MODEL
First, we briefly illustrate our simulation model (more details are provided in Appendix A).

A. Outline
We model a PbSe NC using a polybead structure. The beads forming the NC are disposed to reproduce a rhombicuboctahedron of size 6 nm, that is a typical experimental shape [76,92]. We simulate the dynamics of the NCs using the position-Verlet algorithm [94] to compute the bead's motion. Bead-bead pair potentials are used to reproduce NC-NC short-range facet-specific attractive and repulsive interactions. The distance between beads belonging to the same NC is kept fixed by constraint forces [95]. The solvent is treated implicitly by modeling the NC Brownian motion. External potentials are applied to the NC beads to mimic the interface-adsorption forces experienced by NCs at a fluid-fluid interface.

B. NC-NC interactions
The polybead structure representing each NC is formed by 144 beads (shell beads) disposed to represent the NC surface, see Fig. 1, plus seven beads (core beads) used, see later, for the Brownian motion of the NCs. To model van der Waals and electrostatic-chemical atomic interactions between NCs at almost-contact distance, shell beads of different NCs interact with each other by the Lennard-Jones pair potential, truncated and shifted in r ¼ 2.4σ [see Eq. (A3)], where r is the center-to-center bead-bead distance, σ ¼ 1 nm, and ϵ sets the interaction strength for the bead pair. PbSe NCs are, after their synthesis, dispersed in solution, and ligands, typically oleic acid molecules, are chemisorbed at the NC facets. At this stage, ligands largely screen NC-NC attractive interactions and prevent NC assembly. In the experiments of interest [73,[85][86][87][88], the NC self-assembly is activated by chemically inducing the detachment of ligands from the PbSe NC f100g facets, i.e., where the bonding energy between ligand molecules and NC surface atoms is lowest [73,92]. Hence, during the self-assembly, NCs attract (and attach to) each other only by f100g facets, while NC-NC bonds by f110g and f111g facets are prevented by ligands. To model these NC-NC facet-specific repulsive steric interactions due to ligands, all pairs of shell beads of different NCs involving a red or orange colored bead, see Fig. 1, interact by the soft repulsive pair potential truncated in r ¼ σ R [see Eq. (A5)], with ϵ R ¼ 2 × 10 −22 J, σ R ¼ 2.5 nm when both beads are red or orange, and σ R ¼ 1.25 nm otherwise. These values are order-of-magnitude estimations to reproduce typical soft repulsive potentials between ligand-capped NCs in solvent [96]. Since the NC f100g facets are ligand-free, NCs can attach to each other by opposite f100g facets, performing oriented attachment. The maximum bonding energy E of two NCs attached to each other by f100g facets is expected to be ∼10 −19 J [97]. This value is due not only to NC-NC van der Waals interactions, but also (and mostly) to electrostatic-chemical atomic interactions between surface atoms of the PbSe f100g opposite facets. Note that E is an input parameter of our model and its exact value can be tuned by regulating the value of ϵ in U LJ [Eq. (1)] for the various bead pairs; see Appendix A.

C. NC-solvent interactions
The seven core beads of each NC are placed one in the NC center of mass and the other six at a distance r B from it and in the six h100i directions. To model the NC Brownian motion, the Verlet-type algorithm in Ref. [98] is applied, for each NC, to its central core bead and to three of the remaining six core beads (nonaligned with the central core bead). The parameters of this Verlet-type algorithm and r B are tuned to induce the expected translational and rotational average diffusion of a NC immersed in a typical solvent at constant temperature (room temperature); see Appendix A.

D. Fluid-fluid interface forces
A NC in contact with a fluid-fluid interface experiences a potential energy, depending on its position and orientation, due to the NC interactions with the molecules of the two fluids. Such a potential binds the NC at a certain height at the fluid-fluid interface and drives it to energetically favorable orientations. The specific potentials arising for PbSe NCs at a typical fluid-fluid interface (toluene-air) are predicted in the next section. To reproduce in our coarsegrained MD model a fluid-fluid interface, here identified by the plane z ¼ 0 (where a Cartesian coordinate system x, y, z is introduced), external potentials are applied to the NC beads. To mimic the force bonding a NC to the interface, the potential is applied to the central core bead of each NC, with z c the central core bead z coordinate, u z setting the strength of the force bonding the NC to the interface, z 0 the minimum-U z position of the NC center of mass, z d the maximum distance at which the NC feels the interface (so z d is, roughly, the NC size), and E s , E a , E 1 constants representing the interface-adsorption energy E (see Sec. III) of a single NC, respectively, fully immersed in the fluid below z ¼ 0, fully immersed in the fluid above z ¼ 0, and adsorbed at the interface in its minimum-E position and orientation. To reproduce the forces due to the fluid-fluid interface that drive an adsorbed NC in a specific orientation, we define [see Fig. 2(a)] φ ∈ ½0; π as the polar angle of a NC vertical axis (given by one of its six h100i directions) with respect to the plane z ¼ 0, and ψ ∈ ½0; 2πÞ as the NC internal Euler angle around its vertical axis. Force couples are applied to the NC to induce torques with, respectively, direction orthogonal to the φ, ψ rotational plane, modulus jdU φ =dφj, jdU ψ =dψj, and sign depending on the sign of dU φ =dφ, dU ψ =dψ, where The parameter u p sets the orientation force's strength, and, for simplicity, is assumed constant (although, in principle, it could depend on φ, ψ and be different in U ψ and U φ ). The function ξðz c Þ, which is 1 for z c ¼ z 0 , 0 for jz c j > z d , and a linear interpolation of these values for jz c j < z d , is used to switch on or off the orientation forces when a NC adsorbs or desorbs from the interface. The parameters φ 0 , ψ 0 set the values of φ, ψ with minimum U φ , U ψ . The forces due to U φ , U ψ drive the NC toward the orientation with φ ¼ φ 0 , ψ ¼ ψ 0 . In all the simulations reported in this work, φ 0 , ψ 0 correspond to the NC orientation with a f111g facet parallel to the interface, that is the expected minimum-E orientation for isolated PbSe NCs at a toluene-air interface; see Sec. III. For the rotational symmetry of the NC shape in φ, ψ, different values of φ 0 , ψ 0 correspond to the same NC orientation. For example, there are eight different combinations of φ 0 , ψ 0 corresponding to the orientation of a NC with a specific f111g facet parallel to z ¼ 0 and pointing toward up, one combination for each f111g facet of the NC. In the simulations, the φ 0 , ψ 0 corresponding to the NC minimum-U φ , minimum-U ψ orientation closest to the current NC orientation is chosen, respectively. Cubic and cantellated-rhombicuboctahedron particles, when adsorbed with a f111g facet parallel to the interface, generate hexapolar capillary deformations in the interface height profile inducing capillary interactions between the particles [99,102]. This effect is, however, negligible for the experiments considered here, as we show and discuss extensively in Appendix B 3. Nonetheless, for completeness and future applications, a NC-NC pair capillary potential U c is introduced in our coarse-grained MD model, see Appendix A 4 (although it is irrelevant for the results presented in this work).

E. Thermodynamics considerations
In our coarse-grained MD model, the Helmholtz free energy F of N NCs at a fluid-fluid interface is where T is the (constant) temperature, S N is the entropy of the N NCs, and U tot ≡ U kin þ U pot is the total internal energy of the system, with U kin the total kinetic energy of the NCs and the total potential energy of the NCs. In the definition of U pot [Eq. (7)], it is implicit that U LJ , U R are summed over all bead-bead pairs, U z , U φ , U ψ over all NCs, and U c over all NC-NC pairs. The terms U kin and TS N are included in our model by the NC Brownian motion plus the bead-bead hard interactions (that set the NC shape). The terms U LJ , U R represent the contribution to U pot due to the NC-NC interactions. The interactions of the NCs with the fluid molecules (and the interactions of the fluid molecules between themselves) are implicitly included in our model by implementing the Brownian motion of the NCs. The external potentials U z , U φ , U ψ , U c represent the effects due to the different interactions that NCs at fluid-fluid interfaces have with the molecules of the two different fluids. We call interface-adsorption energy the total energy E ≡ U z þ U φ þ U ψ þ U c due to these "asymmetric" interactions that NCs, when adsorbed at fluid-fluid interfaces, experience with the different molecules of the two different fluids forming the interface (in this definition of E, it is implicit that U z , U φ , U ψ are summed over all NCs, U c over all NC-NC pairs). In Sec. III (and in Appendix B) we (a) Sketch of the configuration z c , φ, ψ of a PbSe NC with respect to the toluene-air interface plane z ¼ 0, where z c is the z coordinate of the NC center of mass, φ ∈ ½0; π is the polar angle of the NC vertical axis v φ (corresponding to one of the NC h100i directions), and ψ ∈ ½0; 2πÞ is the NC internal Euler angle around v φ (and is defined by the vectors v ψ and w, see Appendix A 3). (b)-(d) Interface-adsorption energy E [Eqs. (8)-(10), (B1)] of a PbSe NC, modeled as a rhombicuboctahedron of size 6 nm, at a toluene-air interface, with respect to φ, ψ, z c , computed using the numerical approach in Refs. [99][100][101], with E ¼ 0 corresponding to the NC desorbed from the interface and immersed in toluene. The surface tensions, with toluene and air, assigned to the various NC facets take into account that ligands (oleic acid molecules) are chemisorbed at the NC f110g, f111g facets, while f100g facets are ligand-free (see text). (b) Eðz c Þ for a NC oriented with, respectively, a f100g and a f111g facet parallel to z ¼ 0.
(c),(d) EðφÞ and Eðψ) plotted, respectively, for various values of ψ and φ, and minimized with respect to z c . For the rotational symmetry of the NC shape in φ and ψ, E is shown only for φ ∈ ½0; π=2, ψ ∈ ½0; π=4. The black curves in (b)-(d) show, for various values of u z , u p , the potentials U z , U φ , U ψ [Eqs.
(3)-(5), here opportunely shifted] used in our coarsegrained MD model to represent the interface-adsorption energy E of a PbSe NC at the toluene-air interface. (e),(f) 3D view, close to the NC, of the toluene-air interface minimum-E shape (blue grid, with the toluene the fluid below), computed by the numerical method of Soligno et al. [99][100][101], for the metastable and minimum-E orientation of the NC, i.e., with a f100g and a f111g facet parallel to z ¼ 0, respectively. discuss in more detail the thermodynamic origin of E, and we predict E for the experimental systems of interest (i.e., PbSe NCs at toluene-air interfaces). Since our simulations are performed at constant T, N, and system volume, the NCs will evolve toward a configuration with lower F [Eq. (6)], until a stable minimum in F is reached. A reduction of F (i.e., an energy bonus) can be achieved by increasing S N and/or by decreasing U tot . When the NCs self-assemble, they are evolving toward more ordered structures, so they are decreasing S N (except at very high packing fraction [103]), which corresponds to an increment in F (i.e., an energy penalty). Hence, the NCs self-assemble only if U tot decreases (more than the increment in −TS N ).

III. INTERFACE-ADSORPTION POTENTIALS FOR PbSe NCs
In this section, we compute the potential energy of an isolated PbSe NC at the typical interface toluene-air with respect to the NC configuration z c , φ, ψ at the interface plane, see Fig. 2(a), where a Cartesian coordinate system x, y, z is introduced such that the flat toluene-air interface coincides with the plane z ¼ 0, the toluene is in the semiplane z < 0 and air in the semiplane z > 0. The ultimate purpose of the calculations reported in this section is to estimate the interface-adsorption parameters in our coarse-grained MD model [precisely, in Eqs. (3)-(5)].

A. Macroscopic sharp-interface model
We use a macroscopic model, treating toluene and air as homogeneous fluids, the toluene-air interface as a possibly curved 2D surface with no thickness and equilibrium height profile z ¼ hðx; yÞ, with hðx; yÞ ¼ 0 when no NC is adsorbed, and the PbSe NC as a homogeneous solid with rhombicuboctahedron shape (with two opposite f100g facets at distance 6 nm). Diffuse-interface effects, not included in this approach, become relevant only for particle sizes below 3-4 nm [104]. We assume constant the temperature T, the chemical potentials μ t , μ a of toluene and air molecules, respectively, and the total volume V ¼ V t þ V a þ V NC of the system, with V a , V t the (nonconstant) volumes of air and toluene, respectively, and V NC the (constant) volume of the NC. The convenient thermodynamic potential E for this toluene-air-NC system, henceforth called interface-adsorption energy, is with U the total internal energy of the system, and S t , N t , S a , N a the entropy and total number of molecules of toluene and air, respectively. Gravity is not included in Eq. (8), since it is negligible for the systems of interest. With respect to the toluene and air molecules, E is a grand potential (also called Landau potential). For extensiveness, the grand potential of a homogeneous fluid can be written as −PV, with P the fluid pressure and V the fluid volume.
Hence, E can be written as with P a , P t the pressure of air and toluene, respectively, and U i corresponding to the remaining terms of the internal energy U, due to the toluene-air, NC-toluene, and NC-air interfaces (i.e., the terms not included in the grand potentials of toluene and air molecules). The interfaceadsorption energy E, for a fixed position of the NC, is minimum at equilibrium, since T, V, μ t , μ a are constant. As proven in Ref. [100], the toluene-air interface shape hðx; yÞ that minimizes E, hence the equilibrium hðx; yÞ, is the solution of the Young-Laplace equation, with Young's law for the contact angle holding along the three-phase tolueneair-NC contact lines. With respect to the NC, the potential E is an internal energy. In our coarse-grained MD model, E is represented by the potentials U z , U φ , U ψ , defined in Eqs.
(3)-(5). In the remainder of this section, we compute E for a single-adsorbed PbSe NC at the toluene-air interface, and from it we estimate the parameters for the potentials U z , U φ , U ψ . Since the total volume V is fixed, −P a V a − P t V t in Eq. (9) can be rewritten as −ΔPV t (plus a constant), where ΔP ≡ P t − P a . In our model, we assume a flat toluene-air interface far away from the NC, implying that ΔP ¼ 0 [101], so E is equal to the internal energy U i due to the toluene-air, NC-toluene, and NC-air interfaces. Writing U i explicitly for a single PbSe NC at the toluene-air interface, E is [100] E ¼ γS þ To different facets of the NC we assign different surface tensions with air and toluene. To all the NC f111g and f110g facets, assumed covered by ligands (typically, oleic acid molecules), we assign surface tension 0.018 N=m with air, corresponding to the hexane-air surface tension, and zero with toluene, since ligands are soluble in toluene. Hence, for the f111g and f110g facets, cos θ ≡ ½γ ðaÞ k − γ ðtÞ k =γ ≈ 0.64, where θ is Young's contact angle as defined by Young's law [101]. To all the NC f100g facets, assumed instead ligands-free, we assign cos θ ¼ 0.3, to account for the energetically less-favorable interactions of these facets UNDERSTANDING THE FORMATION OF PbSe … Phys. Rev. X 9, 021015 (2019) 021015-5 with toluene (compared to ligand-toluene interactions). In Eq. (10), the values of S, W ðaÞ k , W ðtÞ k (k ¼ 1; …26) depend on the NC configuration z c , φ, ψ at the interface, and on the interface equilibrium shape h (which also depends on z c , φ, ψ). We use the numerical approach introduced by Soligno et al. [99][100][101], where the equilibrium interface height profile hðx; yÞ is computed with respect to z c , φ, ψ, and then E is obtained. More details on this method are reported in Appendix B 1. Note that, as a first-order approximation, Eðz c ; φ; ψÞ could be computed assuming hðx; yÞ flat, i.e., neglecting capillary deformation effects. However, as shown in Ref. [102], this can lead not only to overestimations of the energy, but even to erroneous equilibrium orientations of the NC.
In Figs Fig. 2(b), Eðz c Þ is plotted for a few fixed orientations (i.e., fixed φ; ψ) of the NC and compared with the potential U z [Eq. (3)], used in our MD model to approximately represent Eðz c Þ, i.e., to induce the forces that keep a NC at the interface plane. As shown, the parameter u z , setting the strength of these forces, is expected to be ≃2 × 10 −20 J for a NC oriented with a f111g facet parallel to z ¼ 0 and ≃5 × 10 −20 J for a NC oriented with a f100g facet parallel to z ¼ 0. In our MD model, for simplicity, we keep u z constant with respect to the NC orientation. In Figs. 2(c) and 2(d), the energy EðφÞ and EðψÞ, minimized over z c , is plotted for various fixed ψ and φ, respectively, and compared with the potentials U φ , U ψ [Eqs. (4) and (5)], used in our MD model to approximately represent Eðφ; ψÞ, i.e., to induce the forces that drive an adsorbed NC toward a certain orientation. As shown, the parameter u p , setting the strength of these forces, is expected to be approximately between 2 × 10 −19 J and 4 × 10 −19 . Interestingly, the results with our MD model (see Sec. IV) show that the NC self-assembly can be directed toward either a square or a silicene-honeycomb superstructure solely by tuning u p (with realistic values for all the other parameters of the MD model). In particular, the turning value for u p is ≃2.5 × 10 −19 (see Sec. IV), i.e., close to its expected value. In the experiments, the self-assembly of PbSe NCs often results in coexistence between square and silicene-honeycomb superstructures, corroborating our predictions. As shown in Fig. 2(f), the NC orientation with minimum E is with a f111g facet parallel to the interface plane z ¼ 0. The energy Eðφ; ψÞ indicates also a metastable orientation of the NC with a f100g facet parallel to z ¼ 0; see Figs. 2(c) and 2(e). Since the energy barriers are small (∼10 −20 J), in our MD simulations we neglect this metastable orientation, defining the interface-adsorption potentials U φ ; U ψ to be minimum only for the NC orientation with a f111g facet parallel to z ¼ 0. Despite this, we show that in certain cases, see Sec. IV, the NCs prefer to orient with a f100g facet parallel to z ¼ 0, since the energy penalty in the free energy F [Eq. (6)] due to the orientational potentials U φ , U ψ is compensated by an energy bonus (greater in absolute value) in the electrostatic-chemical bonds between NCs attached by f100g facets, i.e., in U LJ .

C. PbSe NC fully covered by ligands
In Appendix B 2, we show analogous calculations to those presented in Fig. 2, but for cos θ ¼ 0.64 on all NC facets, corresponding to a NC fully covered by ligands. In this case, multiple metastable orientations of the NC with energy barriers ∼10 −20 J are found, suggesting that, when fully covered by ligands, NCs do adsorb at the interface, but, essentially, with random orientation. This is confirmed by experimental observations [87], where a monolayer of randomly oriented PbSe NCs was observed. Then, in the same experiments [87], a closer-packed monolayer of NCs oriented with a f100g facet parallel to the interface was observed at a later stage. Our interpretation is that ligands were still chemisorbed at (mostly) all the PbSe NCs facets when the monolayer of randomly oriented NCs was observed, while they were (mostly) detached from the NCs f100g facets when the NCs were found oriented. As shown in this section, the expected equilibrium orientation, at the toluene-air interface, of isolated PbSe NCs covered by ligands only on their f110g, f111g facets is with a f111g facet parallel to the interface plane. However, short-range electrostatic-chemical attractive interactions between ligands-free f100g facets of closeby NCs can affect this orientation (see the results with our MD model in Sec. IV), forcing the NCs to have a f100g facet parallel to the interface and resulting in the NCs square superstructures observed in Ref. [87].

D. Considerations on many-particle effects
The free energy F [see Eq. (6), with E represented by U z þ U φ þ U ψ ] can induce a nonminimum-E orientation of the NCs, when many NCs are adsorbed at the interface, because of short-range NC-NC attractive and repulsive chemical interactions (represented by U LJ þ U R ); see Sec. IV. However, the interface-adsorption energy E is expected to be negligibly affected, with respect to the position and orientation z c , φ, ψ of each NC, by multiparticle effects, for the experimental conditions of interest; see Ref. [99]. That is, each adsorbed PbSe NC experiences the same potential Eðz c ; φ; ψÞ we computed for an isolated PbSe NC. So, in our MD model, we are justified in applying the potentials U z , U φ , U ψ , here estimated from the interface-adsorption energy E of a single NC, independently to each interface-adsorbed NC. Capillary interactions between adsorbed NCs, due to the variations of the interface-adsorption energy E with respect to NC-NC reciprocal distances and azimuthal orientations in the interface plane, are also negligible for PbSe NCs adsorbed at a toluene-air interface, as shown and discussed extensively in Appendix B 3.

IV. RESULTS FROM OUR COARSE-GRAINED MD MODEL
In this section, we present the results from our coarsegrained MD model for the self-assembly of PbSe nanocrystals (NCs) of size 6 nm at a fluid-fluid interface.
A. Self-assembly in square and silicene-honeycomb superstructures In Figs. 3 and 4, we present two simulations, with our MD model, of N ¼ 625 NCs, where the NCs self-assemble into, respectively, silicene-honeycomb superstructures (Fig. 3) and square superstructures (Fig. 4). . This value takes into account not only NC-NC van der Waals interactions, but also (and mostly) electrostatic-chemical atomic interactions between PbSe f100g facets at (almost) contact distance, expected since the f100g facets are ligand-free. The parameter u z , setting the strength of the forces due to the interface-adsorption potential U z [Eq. (3)], which are parallel to the z direction and keep the NCs at the interface plane, is u z ¼ 2 × 10 −20 J, i.e., close to the values estimated in Sec. III for PbSe NCs at a toluene-air interface. The NCs' initial configuration (t ¼ 0) is a hexagonal monolayer, with lattice spacing 11 nm, of randomly oriented NCs at the interface.

Initial formation stage
In both simulations (Figs. 3 and 4), the initially dispersed NCs, driven by the potentials U φ , U ψ , U z , orient almost immediately (by t ∼ 0.1 ns) with a f111g facet parallel to the interface plane, that is the minimum-U φ , minimum-U ψ orientation, and keep their minimum-U z height z c at the interface plane. At this stage, the NCs are too far apart to feel the attractive interactions between f100g facets, due to the potential U LJ . Hence, the NCs minimize the total internal potential energy U pot of the system, see Eq. (7), by keeping the configuration z c , ψ, φ with minimum U φ , U ψ , U z (see the graphs of z c , ψ, φ relative to the first snapshot of each simulation, i.e., at t ¼ 0.01 μs). Then, at a later stage (t ≳ 0.01 μs), because of Brownian motion, NCs start bumping into each other, feeling the attractive interactions between f100g facets, due to U LJ . However, as long as the NCs keep the orientation with a f111g facet parallel to the interface plane and the same height z c at the interface plane, they cannot attach by f100g facets, for clear geometrical reasons. Note that NC-NC attachment by f111g or f110g facets would be geometrically possible for NCs in this configuration, but it is prevented by ligands chemisorbed at these NC facets [which, in our model, are represented by the soft repulsive potential U R ; see Eq. (2)].

Formation of silicene-honeycomb superstructures
In Fig. 3, the NCs manage to attach to each other by f100g facets by shifting their height z c at the interface, half of the NCs moving slightly above the interface plane and half of them slightly below, that is, by partially desorbing from the interface (see the two peaks that appear in the z c distribution for increasing time). This shift in z c increases the total potential U z , corresponding to an energy penalty. However, by attaching by f100g facets, the NCs decrease the total potential U LJ more than the increment in U z ; hence, the total internal energy U pot actually decreases. This is confirmed by the energy plots shown in the bottom righthand panel of Fig. 3, where the black line is U pot ðtÞ, the light blue line is U LJ ðtÞ summed over all bead-bead pairs, and the green line is U z ðtÞ summed over all NCs. The NCs' orientation at the interface remains with a f111g facet parallel to the interface plane, i.e., the orientation minimizing U φ ; U ψ , as confirmed by the plot of U φ ðtÞ þ U ψ ðtÞ, summed over all NCs, see fuchsia line in the energy graph, which remains essentially constant (see also the distributions of ψ, φ). As clearly shown also by the simulation snapshots, thanks to this mechanism the NCs aggregate into silicene-honeycomb superstructures, equivalent to those Ligands attached to the NC f110g, f111g facets (preventing NC-NC attachment by these facets) are sketched with black line segments. The graphs at the right of each snapshot show the number n of NCs with a certain value of (green) z c , (red) ψ, and (blue) φ, with z c the NC height at the interface plane and ψ, φ the NC orientation; see Fig. 2

(a). Because of the interface-adsorption potentials
, the NCs are initially in their minimum-E z c , ψ, φ configuration (marked by vertical dotted lines in the z c , ψ, φ distribution graphs), corresponding to NCs all in same plane and with a f111g facet parallel to the interface. As long as the NCs are in this configuration, they are geometrically unable to attach by opposite f100g facets. The attractive forces between opposite f100g facets, due to the potential U LJ [Eq. (1)], drive the NCs out of plane-half of the NCs below, and half of them above, see z c distributions-inducing self-assembly in silicene-honeycomb superstructures. The parameter u p , setting the strength of the interface-adsorption forces orienting the NCs with a f111g facet parallel to the interface plane, i.e., due to U φ , U ψ , is u p ¼ 5.5 × 10 −19 J (for comparison, see Fig. 6, where an equivalent simulation is shown, but for a lower value of u p ). Inset (a) shows an enlargement of a zigzag aggregate of NCs oriented with a f111g facet parallel to the interface plane. Inset (b) shows an enlargement of a large portion of the silicene-honeycomb superstructure formed by the NCs. In particular, defects of this superstructure are highlighted: vacancies (fuchsia circle), lattice's misalignment (see lattices highlighted in blue and in green), trilayer silicene-honeycomb (see bottom left-hand corner of the inset). The graph in the bottom right-hand panel shows, with respect to t, the total potential (light blue) U LJ , i.e., inducing attractions between f100g facets, (green) U z , i.e., keeping the NCs at the interface plane, (fuchsia) U φ þ U ψ , i.e., orienting the NCs at the interface with a f111g facet parallel to the interface plane, and (black) the total internal potential energy U pot ; see Eq. (7). An animation of this simulation is available in the Supplemental Material [105].
, the NCs are initially in their minimum-E z c , ψ, φ configuration (marked by vertical dotted lines in the z c , ψ, φ distribution graphs), corresponding to NCs all in the same plane and with a f111g facet parallel to the interface. As long as the NCs are in this configuration, they are geometrically unable to attach by opposite f100g facets. The attractive forces between opposite f100g facets, due to the potential U LJ [Eq. (1)], tilt the NCs toward non-minimum-E orientations-i.e., with a f100g facet parallel to the interface, see ψ, φ distributions-inducing self-assembly in square superstructures. The parameter u p , setting the strength of the interface-adsorption forces orienting the NCs with a f111g facet parallel to the interface plane, i.e., due to U φ , U ψ , is u p ¼ 1.0 × 10 −19 J (for comparison, see Fig. 5, where an equivalent simulation is shown, but for a higher value of u p ). Inset (a) shows an enlargement of a linear aggregate of NCs oriented with a f110g facet parallel to the interface plane. Insets (b) and (c) show an enlargement of a large portion of the square superstructure formed by the NCs. In particular, defects of this superstructure are highlighted: vacancies (fuchsia circle), lattice's misalignment (see lattices highlighted in blue and in green). The graph in the bottom right-hand panel shows, with respect to t, the total potential (light blue) U LJ , i.e., inducing attractions between f100g facets, (green) U z , i.e., keeping the NCs at the interface plane, (fuchsia) U φ þ U ψ , i.e., orienting the NCs at the interface with a f111g facet parallel to the interface plane, and (black) the total internal potential energy U pot ; see Eq. (7). An animation of this simulation is available in the Supplemental Material [105]. experimentally observed, i.e., where all the NCs have a f111g facet parallel to the interface plane and are disposed in a buckled honeycomb lattice. In inset (a) of Fig. 3, we show an enlargement of a typical precursor of the silicenehoneycomb superstructure, that is, a zigzag aggregate of NCs with a f111g facet parallel to the interface (also experimentally observed [73]). In inset (b) of Fig. 3, we show an enlargement of a large portion of a silicenehoneycomb superstructure formed by the NCs, highlighting defects that have formed: vacancies (i.e., a NC missing from the lattice, see fuchsia circles), misalignment between different domains of the superstructure (see the lattices highlighted in light blue and green), and a trilayer silicenehoneycomb (see the bottom left-hand corner of the inset, where the presence of the extra NCs forming the third layer of the silicene-honeycomb lattice is indicated by the light blue hexagonal lattice). All these defects have been observed experimentally [88]. Their presence in our MD simulations is an additional confirmation that our model correctly captures the key features of the PbSe NC selfassembly in silicene-honeycomb superstructures.

Formation of square superstructures
Differently from the mechanism described for Fig. 3, in Fig. 4 the NCs manage to attach to each other by f100g facets by changing their orientation at the interface, orienting with a f100g facet parallel to the interface plane (see the φ; ψ distributions for increasing time). This change in NC orientation increases the total potential U φ þ U ψ , corresponding to an energy penalty. However, by attaching by f100g facets, the NCs decrease the total potential U LJ more than the increment in U φ þ U ψ ; hence, the total internal energy U pot actually decreases. This is confirmed by the energy plots shown in the bottom right-hand panel of Fig. 4, where the black line is U pot ðtÞ, the light blue line is U LJ ðtÞ summed over all bead-bead pairs, and the fuchsia line is U φ ðtÞ þ U ψ ðtÞ summed over all NCs. All the NCs remain at the same height at the interface plane, keeping the z c value with minimum U z , as confirmed by the plot of U z ðtÞ, summed over all NCs, see green line in the energy graph, which remains essentially constant (see also the distributions of z c ). As clearly shown also by the simulation snapshots, thanks to this mechanism the NCs aggregate into square superstructures, equivalent to those experimentally observed, i.e., where all the NCs are aligned in the same plane and have a f100g facet parallel to the interface plane. In inset (a) of Fig. 4, we show an enlargement of a typical precursor of the square superstructure, that is, a linear aggregate of NCs with a f110g facet parallel to the interface (also experimentally observed [73,85]). In insets (b) and (c) of Fig. 4, we show an enlargement of a large portion of square superstructure formed by the NCs, highlighting defects that have formed: vacancies (see fuchsia circles) and misalignment between different domains of the superstructure (see the lattices highlighted in light blue and green).

B. Assembled-phase diagram
In Figs. 3 and 4 we showed that the NC self-assembly can be drastically affected by solely tuning the parameter u p , regulating the strength of the forces orienting the NCs at the interface with a f111g facet parallel to the interface. Reasonably, on the basis of the formation mechanism for silicene-honeycomb and square superstructures that emerged from the simulations of Figs. 3 and 4, the other two key parameters that should affect the NC selfassembly are E, setting the strength of the forces between f100g facets of different NCs, and u z , setting the strength of the forces keeping all the NCs aligned at the same plane at the interface. To confirm this reasoning, in Fig. 5 we report the assembled phase (square superstructure, silicene-honeycomb superstructure, or disordered, i.e., no self-assembly) obtained from a simulation, with our coarse-grained MD model, of 64 NCs for 2 μs, with respect to the values used for E, u z , u p (using, as NC initial configuration, a hexagonal monolayer with lattice spacing 12 nm at the interface plane z ¼ 0). All the other parameters of our MD model are set to the same values used in the simulations of Figs. 3 and 4 (see Appendix C). Periodic boundary conditions are applied in the x, y directions, with period 96 nm.
As shown in Fig. 5, the obtained NC assembled phase nontrivially depends on E, u z , u p . A low u p (corresponding to weak forces orienting the NCs with a f111g facet parallel to the interface plane) favors formation of square superstructures and, vice versa, a high u p (corresponding to strong forces orienting the NCs with a f111g facet parallel to the interface) favors formation of silicene-honeycomb superstructures. However, silicene-honeycomb superstructures are formed only for low values of u z (corresponding to weak forces keeping the NCs at the interface plane, i.e., preventing partial desorption of the NCs from the interface), while square superstructures are essentially unaffected by u z (for the range considered). Finally, tuning the NC-NC attachment energy E (i.e., the strength of the NC-NC attractive forces) affects the turning values for u z and u p at which different superstructures are observed: a higher NC-NC attraction allows formation of silicenehoneycomb superstructures for higher values of u z , but on the other hand, a higher value of u p is required to prevent formation of square superstructure.

C. Outlook for the experimental synthesis
The range of values considered in Fig. 5 for E, u z , u p represents the realistic intervals in which the experiments [73,[85][86][87][88] are expect to fall, on the basis of the calculations presented in Sec. III for the adsorption at toluene-air interfaces of PbSe NCs (regarding the values of u z , u p ), and on the basis of typical bonding energies between orientedattached nanocrystals [97] (regarding the value of E). As shown by Fig. 5, within these ranges of values, all three phases (disordered, silicene-honeycomb, square) are possible for the experiments. That is, the experiments are very close to the three-phase coexistence point, and slight variations in the experimental conditions can drastically affect the PbSe NC self-assembly, resulting in a different assembled structure. This explains the difficulties encountered in the experiments in controlling and reproducing the NC self-assembly.
In this work, we uncovered the formation mechanism of these superstructures, paving the way toward improved synthesis routes. For example, on the basis of our results, one deduces that the NC self-assembly rate in silicenehoneycomb superstructures can be optimized by increasing the parameters E and/or u p in the experiments, and/or by decreasing u z . The parameter u z can be decreased, e.g., by lowering the solvent-air surface tension, or by using smaller NCs. The parameter u p can be regulated by playing with the solvent-air surface tension and with the size, surface properties, and shape of the synthesized NCs-even slightly different degrees of truncation in the NC shape can considerably affect u p -where NC shape and surface properties can also affect the equilibrium orientation of isolated NCs at the interface. The sharp-interface method of Soligno et al. [99][100][101], here applied in Sec. III, is well suited to predict how these NC properties affect u z , u p . The parameter E can be regulated, e.g., by tuning the NC size or by using a different compound. Clearly, the experimental conditions affecting one parameter in the desired manner might affect another parameter in the opposite way, so care should be used in regulating them. We will carry on systematic studies-based on the model and results presented in this paper-for optimizing the synthesis procedures and will illustrate their findings in future interdisciplinary works with simulations and experiments.

V. UNDERSTANDING THE FORMATION MECHANISM
In this work, we uncovered the formation mechanism of square [73,85,87] and silicene-honeycomb superstructures [85,86,88]; see Secs. IVA 3, IVA 4. Our simulations indicate that silicene-honeycomb self-assembly will occur only if, at some point after the NCs adsorb at the interface, there are deep enough wells in the energy due to NC-NC attachment by f100g facets (i.e., E is low enough) and in the interfaceadsorption energy corresponding to the NC orientation with a f111g facet parallel to the interface plane (i.e., in the potential U φ þ U ψ ), to compensate for the partial desorbment of the NCs from the interface (half of the NCs moving slightly above the interface plane, and the other half slightly below) which causes an energy penalty, because the interface-adsorption energy due to U z increases. Conversely, square self-assembly will occur only if, at some point after the NCs adsorb at the interface, there are deep enough wells in the energy due to NC-NC attachment by f100g facets (i.e., E is low enough) and in the interface-adsorption energy keeping the NCs at the same height at the interface plane (i.e., in the potential U z ), to compensate for the change in orientation of the NCs at the interface (from a f111g to a f100g facet parallel to the interface plane) which causes an energy penalty, because the interface-adsorption energy due to U φ þ U ψ increases. Finally, no NC-NC attachment by f100g facets, hence no self-assembly, is expected if, at any stage after the NCs have adsorbed, the energy well in the NC-NC attachment energy by f100g facets is not deep enough to compensate for either a reorientation of the NCs (required to form the square superstructure) or a partial desorbment of the NCs from the interface (required to form the silicene-honeycomb superstructure), both giving an energy penalty, due to an increase in U φ þ U ψ and U z , respectively.
The most essential results emerging from our work and relevant for understanding the formation mechanism of PbSe square and silicene-honeycomb superstructures are schematized in the next few concise points.
(1) PbSe NCs (size ∼6 nm) fully capped by oleic acid molecules adsorb at a toluene-air interface with energy bonding ≃1.3 × 10 −19 J and, essentially, with random orientation; see Appendix B 2. If many of such NCs adsorb at the interface, a hexagonal monolayer (with randomly oriented NCs) is expected in the close-packed limit, due to the approximately hard-sphere interactions between NCs (since NC-NC attractions are largely screened by ligands).
(a) When isolated, they are bonded to the interface by an energy ≃2.3 × 10 −19 J, and are oriented with a f111g facet parallel to the interface plane; see Sec, III B. (b) When close enough to each other, their behavior is determined by the interplay between NC-NC electrostatic-chemical attractive forces (between NC ligand-free f100g facets) and the interfaceadsorption forces keeping the NCs in plane at the interface and orienting the NCs with a f111g facet parallel to the interface plane. Depending on which forces prevail, the following scenarios occur.
(i) Interface-adsorption forces keeping the NCs in plane at the interface and orienting them with a f111g facet parallel to the interface plane prevail over attractive forces between NC f100g facets. Therefore, NCs cannot attach by opposite f100g facets (because of clear geometrical reasons). A monolayer is formed at the interface-with all NCs oriented with a f111g facet parallel to the interface planedisordered at low NC densities and with hexagonal geometry at high packing fractions. (ii) Attractive forces between NC f100g facets and interface-adsorption forces orienting the NCs with a f111g facet parallel to the interface plane prevail over interfaceadsorption forces keeping the NCs in plane at the interface. The NCs partially desorb from the interface plane (half of the NCs slightly moving up and half of them slightly down) to allow attachment by opposite f100g facets while keeping a f111g facet parallel to the interface plane. Therefore, NCs self-assemble in silicene-honeycomb superstructures. (iii) Interface-adsorption forces keeping the NCs in plane at the interface and attractive forces between NC f100g facets prevail over interface-adsorption forces orienting the NCs with a f111g facet parallel to the interface plane. The NCs orient with a f100g rather than f111g facet parallel to the interface plane (the former orientation is energetically less favorable for an isolated NC) to allow attachment by opposite f100g facets while remaining all in plane at the interface. Therefore, NCs selfassemble in square superstructures.

VI. CONCLUSIONS
In conclusion, introducing a new coarse-grained MD model for the self-assembly of NCs at fluid-fluid interfaces (see Sec. II), we study the self-assembly mechanism of PbSe NCs in square and silicene-honeycomb superstructures, experimentally observed [73,[85][86][87][88]. We show that the final structure is determined by an intricate interplay between interface-adsorption forces orienting and keeping the NCs in plane at the fluid-fluid interface where the self-assembly occurs and short-range electrostaticchemical attractive forces between specific facets of the NCs. Hence, in the formation of these 2D superstructures, a fundamental role is played by the NC interactions with the fluid-fluid interface, in contrast with 3D NC superstructures where, typically, the self-assembly is driven solely by NC-NC entropic and (possibly) attractive interactions. The interface-adsorption forces experienced by the NCs, for typical experimental conditions, are estimated using a macroscopic sharp-interface model; see Fig. 2 (Sec. III). Ligand molecules chemisorbed on the NC surface (weakly on the f100g facets, and strongly on the f110g, f111g facets) also play a fundamental role. Indeed, during the formation process of the PbSe superstructures, ligands are expected to detach from the NC f100g facets, allowing NC-NC attractions only by opposite f100g facets, while NC-NC attachment by f110g, f111g facets is prevented by the ligands still chemisorbed on these facets. The basic mechanism of formation for square and silicenehoneycomb superstructures, arising from the different interplay between the various forces in the system, is illustrated by the simulations shown in Figs As a confirmation that our coarse-grained MD model correctly captures the key features of the dynamics of the NC self-assembly, we show that the formation dynamics of the silicene-honeycomb and square superstructures seems to match the experiments; see Figs. 3 and 4. At the early formation stages of the silicene-honeycomb and square superstructures, we recognize their typical precursors, i.e., zigzag and linear aggregates, respectively (both experimentally observed). In addition, the superstructures formed by the NCs present the same kinds of defects observed experimentally. This suggests that our coarse-grained MD model can also be used to study the formation of these defects and their stability with respect to, e.g., the temperature or other parameters.
We remark that our coarse-grained MD model, here applied to PbSe NCs, can easily be adapted to study the self-assembly at fluid-fluid interfaces of NCs with other shape and interactions. The interface-adsorption potentials, here defined to keep the NCs in the same plane and orient them with a f111g facet parallel to the interface, can easily be modified to reproduce different preferred orientations of the NCs, possibly also including metastable orientations. A graphic-user-interface software implementing our coarse-grained MD model and usable to simulate the NC self-assembly in customizable user-defined conditions will be made (freely) available to the community soon [106]. Our coarse-grained MD approach is computationally much cheaper than fully atomistic or lattice Boltzmann methods, and simulations with ∼1000 NCs at a fluid-fluid interface can easily be performed on a single computer.

ACKNOWLEDGMENTS
The authors acknowledge financial support by the Dutch NWO Physics, FOM program nr. 152 "Engineering Dirac physics in semiconductor superlattices", and the ERC Advanced Grant No. 692691 "FIRST STEP".

APPENDIX A: NC POLYBEAD MODEL
Here we provide more technical details on our coarsegrained molecular dynamics model to simulate the selfassembly of nanocrystals at fluid-fluid interfaces.

Bead dynamics and interactions
Each NC is represented by a polybead structure formed by N a ¼ 151 beads. Each ith bead, with i ¼ 1; …; N a N p , and N p the number of NCs, is formally a pointlike particle, with mass m i and position at the time t given by r i ðtÞ≡ ½x i ðtÞ; y i ðtÞ; z i ðtÞ, where a Cartesian coordinate system x, y, z is introduced. The dynamics of the beads is computed by the position Verlet's algorithm [94]. Constraints forces, to keep beads belonging to the same NC at constant reciprocal distances, are calculated and applied using the method of Ciccotti and Ryckaert [95]. Hence, the basic algorithm to compute the position of the beads at the time t þ Δt, given their positions at the time t and at the time t − Δt, is the following.
(i) Compute the unconstrained bead positions, for i ¼ 1; …; N a N p . (ii) Compute the desired constraint force g i ðtÞ to be applied on each ith bead, given r 0 i ðt þ ΔtÞ, for i ¼ 1; …; N a N p , as illustrated in Ref. [95].
for i ¼ 1; …; N a N p . The time step used in our simulations is Δt ¼ 10 −12 s. The force F i ðtÞ is the total force-applied at the time t to the ith bead-due to the various potentials acting on the beads, that is-see later definitions-the bead-bead pair potentials U LJ , U R , and the external potentials U z , U φ , U ψ , U c used to reproduce the forces on the NCs due to the fluidfluid interface.
The N a beads of each NC are given by seven core beads, used to simulate the NC Brownian motion (see Appendix A 2), plus N a − 7 ¼ 144 shell beads, used to reproduce the NC surface. Shell beads of different NCs interact with each other by the Lennard-Jones pair potential U LJ . Given the generic ith and jth shell beads (i; j ¼ 1; …; N a N p , i ≠ j), the potential U LJ between this pair of beads is defined as where r ≡ jr i − r j j is the distance between the beads, σ is set to 1 nm in our simulations, and c is a constant such that U LJ ð2.4σÞ ¼ 0. The forces f ij and f ji , applied on the ith and jth bead, respectively, and due to the pair potential U LJ between these two beads, are zero if r > 2.4σ, and otherwise, with r ij ≡ r i − r j . In addition, some pairs of shell beads interact also by the soft repulsive pair potential U R , defined as where we set ϵ R ¼ 2 × 10 −22 J in our simulations (and σ R is defined later). If U R is active for the ith and jth bead pair, then the forces f ij and f ji , applied on the ith and jth bead, respectively, and due to the pair potential U R between these two beads, are zero if r > σ R , and otherwise. In Fig. 6, we show a plot of U LJ and U R for illustrative parameters. Periodic boundary conditions can be applied in our model; that is, all the beads in the system are periodically UNDERSTANDING THE FORMATION OF PbSe … Phys. Rev. X 9, 021015 (2019) repeated in the x, y, and z directions with the chosen period. If the minimum image criterium is satisfied, that is, the cutoff radius of each bead-bead pair potential is smaller than half of the periodic box side, then all bead-bead pairs with a given bead can be computed considering only the closest copy of every other bead (since all other pairs are certainly at a distance greater than the cutoff radius). In Fig. 7, we illustrate the disposition of the N a beads used to reproduce a NC. The 144 shell beads are used to reproduce the surface of a NC with rhombicuboctahedron shape and size 6 nm. The seven core beads are placed one in the NC center of mass and the other six at distance r B from it and in the six h100i directions of the NC. For each NC, the mass of the central core bead is 2.114 × 10 −22 kg, the mass of each of the remaining six core beads is 10 −25 kg, and the mass of each of the 144 shell beads is 2 × 10 −24 kg, corresponding to a total mass 5 × 10 −22 kg for the NC. Of course, by using a different number of shell beads and/or placing them in different positions, one can easily tune in our model the NC shape and the accuracy in its representation.
In addition, in our simulations, all shell-bead pairs involving a red or orange colored bead in Fig. 7 interact also by the soft repulsive pair potential U R [Eq. (A5)]. For pairs where both beads are red or orange colored, we assign σ R ¼ 2.5 nm, while for pairs where one bead is red or orange colored and the other is not, we assign σ R ¼ 1.25 nm.

NC Brownian motion
The core beads do not interact by pair potentials, but they are exploited for inducing the Brownian motion of the NCs. For this, the Verlet-type algorithm of Ref. [98] is applied to the central core bead and to three of the remaining six core beads, nonaligned with the central core bead; see the four green colored beads in Fig. 7. That is, the unconstrained . Each of the 24 black colored beads is placed with the center in one of the 24 vertexes of a slightly cantellated rhombicuboctahedron with side 2.5 nm for the f100g square facets, and distance 5 nm between two opposite f100g facets. The position of the center of each blue, light blue, and red colored bead is obtained by a linear interpolation of the centers of two black beads belonging to the same facet. The position of the center of each orange colored bead is obtained by a linear interpolation of the centers of the two closest blue colored beads. The seven core beads, exploited for the Brownian motion of the NC and here colored in green or gray, are placed one in the NC center of mass and the other six at distance r B from it (in these plots, r B ¼ 3 nm) and in the six h100i directions of the NC, respectively. Note that, since the core beads do not interact by bead-bead pair potentials, their position does not affect the shape of the NC even if they are placed outside the NC surface formed by the shell beads. The dynamics of the four green colored core beads follows the Verlet-type Brownian algorithm in Eq. (A8), while the dynamics of the three gray colored core beads, like for all the shell beads, follows the standard position-Verlet algorithm in Eq. (A1). position r 0 i of these four beads, with i the index referring to the bead, is not computed as described in step (i) of the algorithm in [Eq. (A1)], but as follows: for ω ¼ 1, 2, 3, where k B is the Boltzmann constant, T is the temperature expressed in kelvin (in our simulations we set T to room temperature), and ξ iω , for any given time, is an independent random number obtained from a Gaussian distribution with zero mean and unitary variance. The values of the (friction) parameter λ in Eq. (A8) and of r B (i.e., the distance of the side core beads from the central core bead) are tuned to induce the expected translational and rotational average diffusion of a NC immersed in a typical solvent at room temperature. Assuming that a NC approximately behaves like a sphere of diameter D, its translational and rotational diffusion coefficients are [107] respectively, where η is the solvent viscosity. The mean squared displacement (in 3D) of a NC is expected to fulfill the Stokes-Einstein equation [107], i.e., with rðtÞ the NC center-of-mass position at the time t.
The mean squared angular displacement (in 3D) of a NC is expected to fulfill the Stokes-Einstein-Debye relation [107], i.e., with φ the angular displacement of a given direction of the NC after a time t. In our simulations, we set λ ¼ 3 × 10 −11 kg s −1 in Eq. (A8), and we place the side core beads at a distance r B ¼ 3.3 nm from the central core bead. Using these values (heuristically found), the obtained translational and rotational diffusion of the NC fairly matches, see Fig. 8, the diffusion predicted by Eqs. (A14) and (A15), at room temperature, for D ¼ 6 nm, and for η ¼ 5.6 × 10 −4 kg s −1 m −1 , that is, the viscosity of toluene at room temperature. Note that if r B is such that the core beads are placed outside the NC surface formed by the shell beads, the actual shape of the NC is not affected, since the core beads do not interact by bead-bead pair potentials. The remaining three core beads (gray colored beads in Fig. 7) follow the non-Brownian dynamics in Eq. (A1), and are introduced in the model to ensure a distribution of the mass with cubic symmetry with respect to the NC central core bead. Thus far we have assumed the NCs as immersed in a homogeneous fluid. In the next sections, a fluid-fluid interface is introduced in our MD model, by defining external potentials for the NC beads that mimic the forces, due to the interface, experienced by adsorbed NCs. Regarding the NC Brownian motion, we assume that a NC has approximately the same rotational and translational diffusion coefficients when at the fluid-fluid interface. Hence, we keep the same procedure and parameters FIG. 8. Using our coarse-grained MD model (see text), we consider here a single NC in bulk (i.e., no fluid-fluid interface is introduced) and measure the NC average translational and rotational diffusion with respect to the simulation time t, for r B ¼ 3.3 nm (see Fig. 7) and λ ¼ 3 × 10 −11 kg s −1 in Eq. (A8). In (a), the red line shows the squared displacement of the NC center of mass with respect to t, averaged over 1000 simulations. The blue line is the analytical prediction [Eq. (A14)] for a sphere of diameter 6 nm immersed in toluene at room temperature. In (b), the red line shows hφ 2 ðtÞi, that is, φ 2 [with φ the polar angle of the NC vertical axis in the plane z ¼ 0, see Fig. 2(a)], with respect to t and averaged over 1000 simulations, setting φ ¼ 0 as initial configuration. After a time long enough, hφ 2 ðtÞi becomes constant, since, as expected, hφðtÞi (light blue line) goes to π=2. Before this regime takes place, hφ 2 ðtÞi is expected to be linear in t, as predicted by Eq. (A15). The black line is the analytical prediction [Eq. (A15)] for a sphere of diameter 6 nm immersed in toluene at room temperature. The inset is an enlargement of the plot in the initial linear regime of hφ 2 ðtÞi, i.e., for t ≲ 0.05 μs in our case. illustrated in this section to induce the NC Brownian motion when the NCs are at a fluid-fluid interface. The interface effects on the NC dynamics are already taken into account by introducing the external potentials on the beads. Finally, note that hydrodynamics interactions between NCs are not included in our approach, since these are expected to be negligible in the experiments of interest.

Single-NC interface-adsorption potentials
The remaining fundamental ingredient to correctly reproduce in our coarse-grained MD model the selfassembly of PbSe NCs in quasi-2D superstructures is the fluid-fluid interface where the NCs are adsorbed (typically, the interface between the solvent where the NCs are initially dispersed and air). We apply external potentials to the NC beads to mimic the forces due to the fluid-fluid interface.
We assume that the plane z ¼ 0 is the fluid-fluid interface, where a 3D Cartesian coordinate system x, y, z is introduced, with the semiplane z < 0 being one fluid and the semiplane z > 0 being the other fluid. For simplicity, we consider here a single NC. When more NCs are adsorbed at the interface, the same definitions can be applied independently for each NC (Ref. [99] shows that, for the experiments of interest, single-NC interfaceadsorption properties are negligibly affected by multiparticle effects). Given the N a beads forming the NC, labeled by i ¼ 1; …; N a and having position r i ¼ ðx i ; y i ; z i Þ, we indicate by i ¼ 1 the core bead in the NC center of mass (central core bead), and by i ¼ 2; …; 7 the remaining six core beads of the NC (side core beads). The i ¼ 2 and i ¼ 5 side core beads are aligned with the central core bead, and analogously the i ¼ 3 and i ¼ 6, and the i ¼ 4 and i ¼ 7. The three lines passing through these three pairs of side core beads are orthogonal to each other, see Fig. 7, and represent the six h100i directions of the NC. All the other NC beads, i.e., with i ¼ 8; …; N a , are the shell beads. We define external potentials only for the i ¼ 1, i ¼ 2, i ¼ 3, i ¼ 5, and i ¼ 6 beads. The forces exerted on these beads, thanks to the constrained dynamics, are transmitted to the whole NC polybead structure.
To the central core bead of the NC, i.e., with position ðx 1 ; y 1 ; z 1 Þ, we apply the potential U z defined in Eq. (3).
The potential U z represents the interface-adsorption forces binding a NC to the interface, i.e., keeping the NC close to the interface, with z 0 the NC U z height at the interface and u z setting the force strength. The parameter z d represents the maximum distance at which the NC feels the interface (so z d is roughly the NC size). For simplicity, we assumed that U z depends only on the NC center-of-mass distance from the interface, and we used a simple parabolic equation to represent this dependence. More complex and accurate expressions, which, e.g., are not symmetric with respect to z c ¼ z 0 and/or also include the NC orientation, can easily be introduced in our model. To the side core beads of the NC we apply the external potentials inducing the interface-adsorption forces driving the NC toward a specific orientation at the interface. During each time step of a simulation with our MD model, first the orientation of the NC with respect to the interface plane (i.e., z ¼ 0) is computed. Then, given the NC orientation, the interface-adsorption potentials U φ and U ψ are applied to the various side core beads of the NC. To define the orientation of a NC, with respect to the interface plane, we use the angles φ ∈ ½0; π and ψ ∈ ½0; 2πÞ, i.e., the polar angle between the NC vertical axis and the interface plane, and the internal Euler angle of the NC around its vertical axis, respectively; see Fig. 2(a). The vertical axis of the NC is defined as the vector v φ going from the central core bead to the side core bead with position ðx 2 ; y 2 ; z 2 Þ; that is, Then, the angle φ of the NC is Given the side core bead with position in ðx 3 ; y 3 ; z 3 Þ, which is not aligned to ðx 1 ; y 1 ; z 1 Þ and ðx 2 ; y 2 ; z 2 Þ, we define v ψ ≡ ðx 3 ; y 3 ; z 3 Þ − ðx 1 ; y 1 ; z 1 Þ; ðA21Þ and w ≡ ðw x ; w y ; w z Þ as the vector orthogonal to the NC vertical axis v φ and belonging to the plane identified by the vectors (0,0,1) and v φ . Therefore, w · v φ ¼ 0 and w · ½ð0; 0; 1Þ × v φ ¼ 0. Adding the condition jwj ¼ 1, we obtain where we choose the positive sign for w z (and so the negative sign for w x and w y ). Given the angle ψ Ã between w and v ψ , that is, the angle ψ of the NC is defined as Note that, if v φ is parallel to (0,0,1), then w, and so ψ, is not defined. This is not a problem, since, for φ ¼ 0, the NC has a f100g facet parallel to the interface plane z ¼ 0, so its orientation with respect to the interface plane is the same for any ψ; see Fig. 2(a). Once the values of φ and ψ for the NC are computed at a given time step, the external potentials U φ and U ψ are applied to the NC side core beads to induce rotations (by applying force couples) in the rotational planes of φ and ψ. This drives the NC toward orientations with minimum U φ and minimum U ψ . To both the side core beads with position ðx 2 ; y 2 ; z 2 Þ and ðx 5 ; y 5 ; z 5 Þ, i.e., aligned with the central core bead along the direction of NC vertical axis v φ , we apply the rotational potential, with u p setting the force strength, φ 0 the minimum-U φ value of φ for the NC, and ξðz c Þ defined as and representing the dependence of the interfaceadsorption forces due to U φ on the NC distance from the interface. For simplicity, we assumed ξðz c Þ independent from φ and ψ, and set it as a linear function going from zero, when jz c j ≥ z d , to one, when z c ¼ z 0 . Given the vector v φ ¼ ðv φx ; v φy ; v φz Þ defined in Eq. (A19), the force f ¼ ðf x ; f y ; f z Þ due to the potential U φ =2 applied to the side core bead in ðx 2 ; y 2 ; z 2 Þ is while the force due to the potential U φ =2 applied to the side core bead in ðx 5 ; y 5 ; z 5 Þ is −f . Therefore, we are applying a force couple, i.e., a torque, to the NC which, in total, induces a rotation of the NC around its center of mass and in the φ rotational plane. The total potential energy stored by the NC because of its orientation φ at the interface is U φ ðφÞ [Eq (4)]. Note that the forces induced by U φ [Eq. (A27)] are singular in φ ¼ 0 and φ ¼ π, i.e., when v φx ¼ v φy ¼ 0. Indeed, in these cases, the plane of rotation of the angle φ is undefined, since v φ is aligned to z, so the plane passing through v φ and z is undefined. Therefore, these cases need to be treated separately. However, since we consider a potential U φ ðφÞ which is maximum in φ ¼ 0 and φ ¼ π, we simply set the forces to zero. To both the side core beads with position ðx 3 ; y 3 ; z 3 Þ and ðx 6 ; y 6 ; z 6 Þ, i.e., aligned with the central core bead but not along the direction of the NC vertical axis v φ , we apply the rotational potential, with u p setting the force strength, ψ 0 the minimum-U ψ value of ψ for the NC, and ξðz c Þ defined in Eq. (A26) and representing the dependence of the interface-adsorption forces due to U ψ on the NC distance from the interface. Therefore, given the vectors v ψ ¼ ðv ψx ; v ψy ; v ψz Þ and w ≡ ðw x ; w y ; w z Þ defined in Eqs. (A21) and (A22), the force f ¼ ðf x ; f y ; f z Þ due to the potential U ψ =2 applied to the side core bead in ðx 3 ; y 3 ; z 3 Þ is where the sign AE is plus if ψ < π and minus otherwise. The force due to the potential U ψ =2 applied to the side core bead in ðx 6 ; y 6 ; z 6 Þ is −f . Therefore, we are applying a force couple, i.e., a torque, to the NC which, in total, induces a rotation of the NC around its center of mass and in the ψ rotational plane. The total potential energy stored by the NC because of its orientation ψ at the interface is U ψ ðψÞ [Eq. (5)]. Note that, if φ ¼ 0 or φ ¼ π, the angle ψ and the vector w are not defined. In this case, the forces induced by U ψ are set to zero, since there cannot be any potential in ψ when φ ¼ 0 or φ ¼ π (as the NC orientation is the same for any ψ when φ ¼ 0 or φ ¼ π). Note also that the forces induced by U ψ [Eq. (A27)] are singular in ψ ¼ 0 and ψ ¼ π, i.e., when w ¼ AEv ψ , because the rotational plane of the angle ψ, i.e., the plane passing through w and v ψ , is undefined. Therefore, these cases need to be treated separately. However, since we consider a potential U ψ ðψÞ which is maximum in ψ ¼ 0, ψ ¼ π, and ψ ¼ 2π, we simply set the forces to zero. In our simulations, the orientation of a NC with minimum U φ ðφÞ and minimum U ψ ðψÞ is with a f111g facet parallel to the interface plane. This is implemented by setting For simplicity, we assumed u p to be the same parameter in U φ ðφÞ [Eq. (A25)] and U ψ ðψÞ [Eq. (A28)], and we assumed U φ ðφÞ and U ψ ðψÞ independent from ψ and φ, respectively, and used a simple parabolic equation to represent them. More complex and accurate expressions, e.g., including a dependence from both φ and ψ, can in principle be introduced in the method. We remark that the forces induced by U φ and U ψ are approximations meant to capture the key features of the forces experienced by a NC at a fluid-fluid interface. Note that, although φ ∈ ½0; π and ψ ∈ ½0; 2πÞ provide a one-to-one map to all the possible orientations of a NC with respect to the interface plane z ¼ 0, see Fig. 2(a), the rotational planes of φ and ψ are not necessarily the same if we consider certain symmetric orientations of the NC (consider, e.g., a NC with a f100g facet parallel to the interface plane). This slight inconsistency occurs because we are using only two planes of rotation (defined by φ and ψ) to change the orientation of a 3D object. To avoid this, one should introduce a third plane of rotation, orthogonal to the rotational planes of φ and ψ, and then define a rotational potential also in this plane. Anyway, for the orientations of a NC with a f111g facet parallel to the interface plane, which is the NC orientation with minimum-U φ and minimum-U ψ considered in our simulations, this problem does not occur.

NC-NC pair capillary potential
As predicted in Refs. [99,102], cubes and cantellated rhombicuboctahedrons, when adsorbed at a fluid-fluid interface with a f111g facet parallel to the interface plane, generate hexapolar capillary deformations in the interface height profile that induce capillary interactions between the particles. We show in Appendix B 3 that this effect is negligible for the self-assembly of PbSe NCs in the experiments of interest. Nonetheless, for completeness and future applications, we introduce in our coarse-grained MD model the NC-NC pair potential U c , with the form of the expected interaction energy between two hexapolar capillary deformations [108], to include NC-NC capillary interactions, at least in an approximate form (i.e., neglecting multiparticle effects and considering only the case of NCs with a f111g facet parallel to the interface).
To introduce the NC-NC capillary pair potential U c , we assume here that only two NCs, each with N a beads, are at the interface. If more NCs are present, the same definitions are applied to each pair of NCs. The index i ¼ 1; …; 2N a is used to label the 2N a beads in the system, each with position r i ¼ ðx i ; y i ; z i Þ, with the first N a beads (i ¼ 1; …; N a ) belonging to the first NC and the remaining N a beads (i ¼ N a þ 1; …; 2N a ) belonging to the second NC. The position of the center of mass, i.e., of the central core bead, of the first and second NC is, respectively, r 1 ¼ ðx 1 ; y 1 ; z 1 Þ and r N a þ1 ¼ ðx N a þ1 ; y N a þ1 ; z N a þ1 Þ. We introduce φ p , ψ p , α p , with p ¼ 1, 2, where φ p and ψ p are the angles φ and ψ, as defined in Eqs. (A20) and (A24), of the first and second NC, respectively, and α 1 , α 2 are the azimuthal orientations in the plane z ¼ 0, and measured from the direction joining the NC centers of mass, of the vertical axis of the first and second NC, respectively. Given the vertical axis of the first NC, i.e., v φ 1 ≡ ðx 2 ; y 2 ; z 2 Þ − ðx 1 ; y 1 ; z 1 Þ; ðA31Þ the vertical axis of the second NC, i.e., v φ 2 ≡ ðx 2þN a ; y 2þN a ; z 2þN a Þ − ðx 1þN a ; y 1þN a ; z 1þN a Þ; ðA32Þ and the vector v D defining the line joining the centers of mass of the two NCs, i.e., we consider their projections v φ1xy ≡ ðv φ1x ; v φ1y ; 0Þ, v φ2xy ≡ ðv φ2x ; v φ2y ; 0Þ, v Dxy ≡ ðv Dx ; v Dy ; 0Þ on the z ¼ 0 plane, i.e., v φ 1xy ¼ ðx 2 − x 1 ; y 2 − y 1 ; 0Þ; Given α Ã 1 and α Ã 2 , defined as we define α 1 and α 2 as The total energy U c of the two NCs, each inducing a hexapolar capillary deformation in the interface, is theoretically expected to be of the form [108] U c ðD; α 1 ; is the distance between the NC centers of mass in the z ¼ 0 plane, where the factor cosð3ΔαÞ takes into account the alignment of the hexapoles (such that the NC-NC force due to U c is repulsive when a capillary rise of the hexapole induced by a NC overlaps with a capillary depression of the hexapole induced by the other NC and attractive when two capillary rises or two capillary depressions overlap), u c is a parameter setting the interaction strength, and the factor H is used to switch on and off, during the simulation, the pair capillary interaction, depending on whether both NCs are in the orientation inducing a hexapolar capillary deformation or not. The phase α φ 1 ;φ 2 is and is introduced because, when φ ¼ 0.7π, all the signs of the capillary deformations of the hexapolar deformation induced by the NC are inverted from to the case φ ¼ 0.3π (i.e., rises and depressions in the interface height profile are inverted). The factor H is defined as where HðxÞ is the Heaviside step function, z 0 is the minimum-U z height of the NC center of mass, φ . These values are estimations to take into account that a NC, when slightly deviating from the orientation with a f111g facet parallel to z ¼ 0, still induces a hexapolar capillary deformation in the interface height profile. Capillary deformations that may arise for NCs in other orientations, i.e., far from the orientation with a f111g facet parallel to z ¼ 0, are neglected since these orientations are nonstable for the experimental systems of interest, as shown in Sec. III.
The capillary forces induced by the pair potential U c ðD; α 1 ; α 2 Þ [Eq. (A41)] on the two NCs have a translational component acting on the distance between the two NCs and a rotational component acting on the angles α 1 and α 2 of the two NCs. The translational component of the force is obtained by applying to both the i ¼ 1 and i ¼ N a þ 1 beads, i.e., the central core beads of the two NCs, the potential U c ðDÞ [Eq. (A41)]. The force f ¼ ðf x ; f y ; f z Þ applied to the i ¼ 1 bead and due to U c ðDÞ is Then, the force applied to the i ¼ N a þ 1 bead and due to U c is −f . To implement the rotational component of the force due to U c , we select two shell beads, for each NC, aligned with the NC center of mass (so that we can apply a force UNDERSTANDING THE FORMATION OF PbSe … Phys. Rev. X 9, 021015 (2019)

021015-19
couple to the NC) and in the same x, y plane, so that the induced rotation occurs in the rotational planes of the angles α 1 and α 2 , i.e., in a plane parallel to z ¼ 0. Given our definition of the shell bead positions for a PbSe NC, see Fig. 7, we can always find two shell beads that fulfill these properties (within 0.02 nm of precision) when the NC has a f111g facet exactly parallel to the interface. When a NC is slightly tilted from this orientation, we select the shell bead pair that fulfills these properties with the best approximation. We now assume the two selected shell beads with such properties are the i ¼ 8 and i ¼ 30 beads for the first NC, and the i ¼ N a þ 8 and i ¼ N a þ 30 beads for the second NC. Their position is, respectively, ðx 8 ; y 8 ; z 1 Þ, ðx 30 ; y 30 ; z 1 Þ, ðx Naþ8 ; y Naþ8 ; z Naþ1 Þ, ðx Naþ30 ; y Naþ30 ; z Naþ1 Þ. To the i ¼ 8 bead we apply a force f ¼ ðf x ; f y ; 0Þ orthogonal to a ≡ ðx 8 − x 1 ; y 8 − y 1 ; 0Þ ≡ ða x ; a y ; 0Þ, and to the i ¼ 30 bead we apply a force −f , such that the total modulus of the torque τ α 1 induced on the first NC by this force couple is This implies that the force f is given by Note that, since the i ¼ 8 and i ¼ 30 beads are aligned with the center of mass of the NC in the same x, y plane, the direction of the torque induced by f and −f on the NC is orthogonal to the plane z ¼ 0. So, the NC rotation occurs, as desired, in the rotational plane of α 1 . The signs AE for f y and ∓ for f x in Eq. (A48), determining the sign of the torque τ α 1 , depend on the sign of ∂U c =∂α 1 : AE is þ (and ∓ is −) if −∂U c =∂α 1 is positive, and AE is − (and ∓ is þ) otherwise. Analogously, for the second NC, we apply to the i ¼ N a þ 8 bead a force f ¼ ðf x ; f y ; 0Þ orthogonal to b ≡ ðx N a þ8 − x N a þ1 ; y N a þ8 − y N a þ1 ; 0Þ ≡ ðb x ; b y ; 0Þ, and to the i ¼ N a þ 30 bead we apply a force −f , such that the total modulus of the torque τ α 2 induced on the second NC by this force couple is also τ α [Eq. (A47)]. Hence, the force f is given by The signs AE for f y and ∓ for f x in Eq. (A49), determining the sign of the torque τ α 2 , depend on the sign of ∂U c =∂α 2 : AE is − (and ∓ is þ) if −∂U c =∂α 2 is positive, and AE is þ (and ∓ is −) otherwise. Note that, since ∂U c =dα 2 ¼ −∂U c =∂α 1 , it follows that τ α 1 ¼ τ α 2 ; hence if τ α 1 acts to increase α 1 , then τ α 2 acts to decrease α 2 , and vice versa.
Using the typical value u c ¼ 5 × 10 −16 J nm 6 (see Appendix B 3), the maximum value for jU c j [Eq. (A41)] with respect to α 1 , α 2 is ≃ − 5 × 10 −23 J at a distance D ≡ r U c ≡ 14.4 nm between the centers of mass of the two interacting NCs. Since this value is definitely negligible compared to the energy scales involved in our MD model, we use r U c as a cutoff radius for the NC-NC pair interaction potential U c [Eq. (A41)]; i.e., we neglect the pair capillary interactions of NC pairs at a distance greater than r U c . When periodic boundary conditions are applied, if r U c is less than half of the periodic box side, then, for each NC, only the NC-NC pairs with the closest copy of every other NC need to be considered.
In all the simulations with our coarse-grained MD model presented in this paper, the pair potential U c is used, setting u c ¼ 4 × 10 −16 J nm 6 . Its effect on the NCs' dynamics is, however, negligible. During each simulation, the total potential U c (i.e., summed over all NC-NC pairs) is essentially always constant, compared to the other potentials U z , U φ , U ψ (summed over all NCs) and U LJ (summed over all bead-bead pairs). Using values of u c a few times greater, and removing short-range attractive forces, we find instead that U c affects the NCs' dynamics, driving the NCs toward the hexagonal-lattice assemblies predicted in Refs. [99,102]. However, for conciseness, we leave this out of this work.

APPENDIX B: NC INTERFACE-ADSORPTION ENERGY CALCULATIONS
Here we report more details on the calculations presented in Sec. III for the adsorption of PbSe NCs at a toluene-air interface and show additional results.

Sharp-interface method
First, we report additional details on the numerical method introduced by Soligno et al. [99][100][101] and used for the interface-adsorption calculations presented in Sec. III and in Appendixes B 2 and B 3.
We consider a macroscopic model where the interface between two homogeneous fluids, here assumed toluene and air, is treated as a 2D possibly curved surface, which we represent by a 2D triangular grid of points. A PbSe NC adsorbed at this interface is modeled as a homogeneous solid with rhombicuboctahedron shape (with distance 6 nm between two opposite f100g facets), and its surface is represented by a 2D closed triangular grid of points. We numerically calculate the energy E [Eq. (10)] of the NCtoluene-air system with respect to the NC position and orientation at the interface. First, the equilibrium shape of the toluene-air interface for the given NC position and orientation is computed, then E is obtained. As a first-order approximation, one could assume the interface height profile always flat and unperturbed by the presence of a particle at the interface, that is the approach first introduced by Pieranski [109]. For example, numerical methods applying this approximation are used in Refs. [26,[110][111][112][113][114][115][116]. However, neglecting capillary deformations can lead to erroneous predictions even for an isolated adsorbed particle; see Ref. [102].
For convenience, the interface-adsorption energy E [Eq. (10)] is shifted by a constant and rewritten in the equivalent form where S and A are the surface areas of the toluene-air interface with and without NC, respectively (so A is a constant), the index k goes over the 26 facets of the NC, and W First, the minimum-E shape of the toluene-air interface grid is computed, given as input a fixed position of the NC surface, defined, see Fig. 2(a), by the position and orientation z c , φ, ψ of the NC with respect to the toluene-air interface plane (corresponding to z ¼ 0, where a Cartesian coordinate system x, y, z is introduced). For this calculation, we use a simulated annealing algorithm (i.e., a Monte Carlo approach) to calculate the position with minimum E [Eq. (B1)] of the points of the toluene-air interface grid. As shown in Ref. [100], the obtained interface shape is the solution of both the Young-Laplace equation and Young's law [101], that is, the equilibrium shape. After the interface minimum-E shape is obtained, S and W ðaÞ k (for k ¼ 1; …; 26) are computed from the position of the interface grid, and, finally, E is obtained from their values; see Eq. (B1). The toluene-air interface area S is obtained by summing the areas of all the triangles forming the air-toluene grid. The area of the portion of the NC kth facet in contact with air, i.e., W ðaÞ k , is obtained by summing the areas of all the triangles of the NC surface grid belonging to the NC kth facet and above the air-toluene grid (with the triangles in between the tolueneair interface split, and only the surface area of their portion above the air-toluene grid included). By repeating this procedure for different values of φ, ψ, z c , the energy landscape Eðφ; ψ; z c Þ is obtained.
The computed minimum-E air-toluene interface shape forms an angle with the NC surface along the three-phase contact line that matches the input parameter θ k (with k index of the facet where the three-phase contact line is considered). The position of the three-phase contact line, for a given NC configuration φ, ψ, z c , is automatically found by minimizing E with respect to the positions of the air-toluene interface grid points; i.e., it is not imposed a priori.
The initial shape of the toluene-air interface grid, i.e., before starting a simulated annealing simulation to minimize E, is the plane z ¼ 0. To simulate a toluene-air interface which is flat far away from the NC, the tolueneair-NC system is enclosed by a wall orthogonal to the plane z ¼ 0 and with Young's contact angle π=2, placed far enough from the NC to avoid finite-size effects. During the simulated annealing, we can either allow volume exchange between air and toluene or not. So, the equilibrium shape of the toluene-air interface is found, respectively, for the minimum-E volume of toluene and air (such that their sum is constant) or with the constraint of fixed toluene and air volumes (defined by the initial interface shape). In Fig. 2(b), where the energy Eðz c Þ is shown for some fixed orientations φ, ψ of the NC, the interface shape minimizing E is found imposing fixed volumes for toluene and air. In Figs. 2(c) and 2(d), where the energy EðφÞ, EðψÞ is shown for some fixed values of ψ; φ, respectively, the toluene-air interface shape minimizing E is computed allowing volume exchange between toluene and air; hence, the obtained energy Eðφ; ψÞ is minimized over z c , i.e., the NC height at the interface plane.

Interface-adsorption energy of PbSe NCs
In Sec. III, we presented results for the interfaceadsorption energy E of an isolated PbSe NC at a toluene-air interface, with respect to the NC position and orientation z c , φ, ψ [see Fig. 2(a)] at the interface plane z ¼ 0. Ligand molecules (typically, oleic acid) were assumed chemisorbed at the NC f111g and f110g facets, while the NC f100g facets were assumed ligands-free. To represent this in the macroscopic model used for calculating E, we assigned a Young's contact angle θ given by cos θ ¼ 0.30 on the NC f100g facets and by cos θ ¼ 0.64 on the NC f110g, f111g facets; see Sec. III.
In this section, we present analogous results, but for a PbSe NC fully covered by ligands; that is, we set cos θ ¼ 0.64 on all the NC facets. In Fig. 9, we show, for such a NC, the energy Eðφ; ψÞ [Eq. (B1)], minimized on z c (i.e., allowing volume exchange between toluene and air when computing the minimum-E shape of the toluene-air interface). In contrast with the energy Eðφ; ψÞ of the NC with f100g facets ligands-free, which presented a deep well (energy barriers ≃4 × 10 −20 J) around the minimum-E orientation (i.e., with a f111g facet parallel to z ¼ 0, see Fig. 2), the energy Eðφ; ψÞ of the NC fully covered by ligands is characterized by several minima with smaller energy barriers (≃2 × 10 −20 J); see Fig. 9. This suggests that, when fully covered by ligands, PbSe NCs at the toluene-air interface have essentially random orientations, since the NC can easily jump from metastable to metastable orientation by thermal motion. Note, however, that PbSe NCs fully covered by ligands are bonded to the air-toluene interface by an energy ∼10 −19 J (see Fig. 9, where the energy level E ¼ 0 corresponds to the NC desorbed from the interface and immersed in toluene); hence, once adsorbed at the interface, they are expected to remain at the interface, even if they can change orientation. For PbSe NCs with f100g facets free from ligands, the energy bonding to the air-toluene interface is even stronger (≃2 × 10 −19 J, see Fig. 2), so, in the experiments, when ligands detach from the f100g facets of the PbSe NCs, the NCs remain strongly bonded to the toluene-air interface.

Capillary interactions of PbSe NCs
We study in this section the capillary interactions between PbSe NCs adsorbed at a toluene-air interface, with the ultimate goal of showing that these interactions are negligible compared to the other forces involved in the NC self-assembly, for the experimental conditions of interest.
In Sec. III and Appendix B 2, we studied the interfaceadsorption energy E [Eqs. (8)-(10), (B1)] for a singleadsorbed PbSe NC at the toluene-air interface, with respect to the NC position and orientation z c , φ, ψ [see Fig. 2(a)] at the interface plane z ¼ 0. If many PbSe NCs are adsorbed at the interface, the energy E with respect to only the position and orientation z c , φ, ψ of each NC is expected to be negligibly affected by multiparticle effects; see Ref. [99]. Hence, the minimum-E values of z c , φ, ψ of each NC, computed for a single-adsorbed PbSe NC, are the same when many PbSe NCs are at the toluene-air interface. However, when many NCs are adsorbed, if they induce capillary deformations in the equilibrium shape of the fluidfluid interface, the interface-adsorption energy E can vary, hence generate interactions, with respect to the reciprocal distances and reciprocal azimuthal orientations in the interface plane between the adsorbed NCs. Indeed, the capillary deformations induced by the adsorbed NCs in the interface shape correspond to a larger surface area of the interface, i.e., the term S in Eqs. (10) and (B1). Therefore, adsorbed NCs inducing capillary deformations can arrange themselves to overlap capillary deformations with the same sign (i.e., either rises or depressions in the interface height profile), decreasing in this way S, and so E. These interactions between adsorbed NCs are called capillary interactions [102].
References [99,102] show that cubes and cantellated rhombicuboctahedrons can adsorb at a fluid-fluid interface with a f111g facet parallel to the interface plane, and in this orientation they induce a hexapolar capillary deformation in the interface height profile. From the capillary interactions induced by such capillary deformations, and taking into account also the entropy of the particles (using approximated analytic expressions), the self-assembly into periodic lattices with honeycomb, hexagonal, and square geometry is predicted. These lattices are different from the NCs' superstructures observed experimentally and predicted in this work with our coarse-grained MD model; for example, the honeycomb lattice is not buckled, the FIG. 9. Interface-adsorption energy E [Eq. (B1)] of a PbSe NC fully covered by ligands, that is with Young's contact angle θ given by cos θ ¼ 0.64 on all the NC facets, with respect to the NC orientation φ, ψ, see Fig. 2 (and minimized on the NC height z c ), at the airtoluene interface plane z ¼ 0, computed using the numerical approach in Refs. [99][100][101]. The level E ¼ 0 corresponds to the NC desorbed from the interface and immersed in toluene. For the rotational symmetry of the NC shape in φ and ψ, E is shown only for φ ∈ ½0; π=2, ψ ∈ ½0; π=4, since outside these intervals is periodically repeated. The insets on the right-hand side show a 3D view, close to the NC, of the toluene-air interface minimum-E shape (blue grid, with the toluene below), computed by the numerical approach of Soligno et al. [99][100][101], for the various metastable-E orientations of the NC, that is, particles in the square lattice are oriented with a f111g, rather than f100g, facet parallel to the interface plane, and in general the particles in these lattices are not attached by f100g facets. As a matter of fact, in Refs. [99,102] only capillary and hard interactions between adsorbed particles at fluid-fluid interfaces are considered, i.e., in the limit of no short-range forces. Short-range facet-specific attractions between NCs, however, are crucial for the self-assembly of PbSe NCs, as shown by the results with our coarse-grained MD model. Given, however, the fair similarity of the honeycomb lattice predicted in Refs. [99,102] with the PbSe NCs silicene-honeycomb superstructure experimentally observed, and that the narrow range of parameters predicted for obtaining the honeycomb lattice was compatible with the PbSe NCs experiments, in Refs. [99,102] it was speculated that NC-NC capillary interactions could play a major role in the self-assembly of NCs at fluidfluid interfaces (although, as stated by Soligno et al., short-range forces needed to be included to verify this hypothesis). We show in this section that NC-NC capillary interactions are actually negligible in the NC self-assembly, at least for the experiments considered here [73,[85][86][87][88]. Note that, in Refs. [99,102], an estimated comparison between the predicted NC-NC capillary interactions and NC-NC van der Waals attractions was provided, showing that, for typical NC sizes, van der Waals attractions are negligible compared to capillary interactions, except at almost-contact NC-NC distance where these two forces are of the same order. However, in the experiments considered here, the f100g facets of PbSe NCs are not covered by ligands; hence, at short range, electrostatic-chemical interactions between surface atoms of opposite f100g facets (∼10 times stronger than van der Waals forces) are also involved. Note also that in Refs. [99,102] it was assigned a Young's contact angle θ given by cos θ ¼ 0 on the whole particle surface, while here we use more realistic surface tensions for the experiments of interest.
We consider in this section PbSe NCs at a toluene-air interface, assigning cos θ ¼ 0.30 on their f100g facets and cos θ ¼ 0.64 on their f110g, f111g facets. As shown in Sec, III, the minimum-E orientation of these NCs is with a f111g facet parallel to the interface plane. When adsorbed in this orientation, each NC induces in the toluene-air interface height profile a hexapolar capillary deformation; see Fig. 10. To estimate the capillary interactions induced by these hexapolar capillary deformations, we consider the hexagonal and honeycomb NCs' periodic lattices predicted in Refs. [99,102], i.e., which minimize E by capillary interactions. To consider a periodic lattice of NCs at the interface, we apply our numerical method for computing the minimum-E interface shape (see Appendix B 1) to a unit cell of the lattice, with periodic boundary conditions applied at the borders of the unit cell to the interface height profile. Each NC in the lattices is oriented with the φ; ψ orientation of minimum E computed for a single-adsorbed NC, i.e., with a f111g facet parallel to the interface plane. The center-of-mass height z c of each NC is aligned on the same plane (parallel to the z ¼ 0 plane), and the minimum-E shape of the toluene-air interface is computed allowing volume exchange between air and toluene. The hexagonallattice unit cell is a hexagon with each side having the same interface height profile of its opposite side. A single NC is in the center of the hexagonal-lattice unit cell. The NC azimuthal orientation in the lattice is shown in the contour plots of Fig. 11. The honeycomb-lattice unit cell is a rectangle with short sides having the same interface height profile, and long sides having the same interface height profiles but shifted of half-side. Each cell of the honeycomb lattice has two NCs, positioned and azimuthally oriented in the lattice as shown in the contour plots of Fig. 11. An exact definition of these two unit cells is reported in Ref. [99].
In Figs. 11(a) and 11(b), we show the pair capillaryinteraction energy E Ã for the honeycomb-and hexagonallattice unit cell, respectively, as computed with our numerical method and with respect to the center-ofmass distance D between two closest-neighbor NCs in the lattice, where E Ã is the interface-adsorption energy E [Eq. (B1)] divided by the number of capillary bonds (between closest neighbors) in the lattice, and shifted to be zero for an infinite lattice spacing. That is, where n is the number of pairs between closest-neighbor NCs in the lattice, N is the number of NCs in the lattice, and  [99][100][101], close to a PbSe rhombicuboctahedron NC (with Young's contact angle θ given by cos θ ¼ 0.64 on the NC f111g, f110g facets and cos θ ¼ 0.3 on the NC f100g facets) in its minimum-E orientation, i.e., with a f111g facet parallel to the interface plane (see Sec. III). In this orientation the NC induces a hexapolar capillary deformation in the interface. On the right, a 3D profile view is shown, with the blue grid being the toluene-air interface (and the toluene below).
which is the expected potential U c between two interacting hexapolar capillary deformations, see Eq. (A41) in Appendix A 4, for cos½3ðα 1 − α 2 þ α φ 1 ;φ 2 Þ ¼ 1. Finally, in Fig. 11(c), we show the energy E Ã [Eq. (B2)] of a honeycomb-lattice unit cell, for a fixed lattice spacing (D ¼ 7.8 nm) and with respect to the azimuthal orientation α in the lattice unit cell of one of the two NCs (keeping fixed the azimuthal orientation in the lattice of the other NC). In Fig. 11(c), we also show a fit of the obtained energy E Ã ðαÞ using gðαÞ ≡ − u c ð7.8 nmÞ 6 cos½3ðα þ π=2Þ; ðB4Þ which is the potential U c [Eq. (A41)] at the NC-NC distance 7.8 nm.
The parameter u c sets the strength of the forces due to U c . The estimated values of u c in the three cases shown in Figs. 11(c)-that is, respectively, u c ¼5.7×10 −16 Jnm 6 , u c ¼ 4.5 × 10 −16 J nm 6 , and u c ¼ 4.6 × 10 −16 J nm 6are, as expected, similar. Note that, in this approach to estimate the pair capillary interactions between PbSe NCs at the toluene-air interface, we are neglecting multiparticle effects. That is, the estimated pair capillary interaction potential U c is actually stronger than in the case of an isolated pair of NCs. Despite this overestimation, the maximum bonding energy between two capillary-interacting NCs, that is, fðDÞ [Eq. (B3)] in the hexagonal lattice at the contact distance D ≃ 6 nm, is ≃1.3 × 10 −20 J, i.e., negligible compared to the total bonding energy E ∼ 10 −19 J expected from the electrostatic-chemical atomic interactions between NCs attached by opposite f100g facets. Hence, capillary interactions are not expected to play a role in the self-assembly of PbSe NCs at the toluene-air interface. Instead, as shown by our coarsegrained MD model, the NC self-assembly is regulated by the interplay between interface-adsorption forces binding and orienting the NCs at the interface and the attractive electrostatic-chemical forces between specific facets of NCs at close distance. The pair capillary potential U c is implemented in our coarse-grained MD model, see Appendix A 4, and used in all the simulations presented in this work setting u c ¼ 4 × 10 −16 J nm 6 . However, U c does not play a role in the NCs' dynamics. Indeed, the total potential U c , i.e., summed over all NC-NC pairs, in all the simulations is essentially always constant with respect to the potentials U z , U φ , U ψ (summed over all NCs) and U LJ (summed over all bead-bead pairs). Even in the simulations where no NC attachment by f100g facets is achieved, see results in Appendix C, the NCs' dynamics are essentially unaffected by U c .

APPENDIX C: ADDITIONAL SIMULATION DETAILS
Here, we report additional simulation details for the results shown in Figs. 3-5 of Sec. IV.
For the interface-adsorption potential U z [Eq. (3)], we set z 0 ¼ −1 nm, z d ¼ 4 nm, E 1 ¼−2 ×10 −19 J, E a ¼ 10 −18 J, E s ¼ 0 (the latest three parameters do not matter for the simulation, they just set the offsets of U z for the energy plots). For the NC-NC pair capillary potential U c introduced in Appendix A 4, see Eq. (A41), we set u c ¼ 4 × 10 −16 J nm 6 , which is close to the value we estimated for PbSe NCs at a toluene-air interface; see Appendix B 3. As extensively discussed in Appendix B 3, for such a value of u c , the NC-NC pair potential U c does not play a role in the NC self-assembly, since the capillary interactions between NCs oriented with , with respect to the center-of-mass distance D between two closestneighbor NCs in a honeycomb-(a) and hexagonal-lattice (b) unit cell (see text), as computed through the numerical approach illustrated in Appendix B, for PbSe NCs of size 6 nm at a toluene-air interface, with Young's contact angle θ given by cos θ ¼ 0.30 on the NC f100g facets, cos θ ¼ 0.64 on the f110g, f111g facets. Each NC is oriented at the interface with a f111g facet parallel to the interface plane, hence generating a hexapolar capillary deformation; see Fig. 10. The black curves show a fit using fðDÞ [Eq. (B3)], which gives u c ¼ 4.5 × 10 −16 J nm 6 for the honeycomb lattice and u c ¼ 5.7 × 10 −16 J nm 6 for the hexagonal lattice. The vertical dotted lines indicate the NC-NC contact distance. (c) Analogous calculation to (b), i.e., for a honeycomb-lattice unit cell, but for a fixed lattice spacing (D ¼ 7.8 nm) and rotating the azimuthal angle α in the lattice of one of the two NCs. The black curve shows a fit using gðαÞ [Eq. (B4)], which gives u c ¼ 4.6 × 10 −16 J nm 6 . The insets show contour plots of the equilibrium interface height profile hðx; yÞ, as computed by the numerical method of Soligno et al. [99][100][101], for some of the unit cells considered.
FIG. 12. For the assembled-phase diagram of Fig. 5, in the case of total bonding energy between two NCs attached by f100g facets given by E ¼ −0.7 × 10 −19 J, we show here a top view of the interface at the end of each simulation (t ¼ 2 μs) for the various values considered of u p , u z (i.e., the parameters setting the strength of the interface-adsorption forces orienting and keeping the NCs at the interface plane, respectively). The black line borders indicate the periodic boundary conditions. The graphs above each snapshot show, for each simulation and with respect to t=μs, the Lennard-Jones potential (light blue line) U LJ =10 −17 J summed over all the bead-bead pairs, the interface-adsorption potentials (green line) U z =10 −17 J and (fuchsia line) ðU φ þ U ψ Þ=10 −17 J, summed over all NCs, and the total potential internal energy (black line) U pot =10 −17 J; see Eq. (7).
FIG. 13. For the assembled-phase diagram of Fig. 5, in the case of total bonding energy between two NCs attached by f100g facets given by E ¼ −1.0 × 10 −19 J, we show here a top view of the interface at the end of each simulation (t ¼ 2 μs) for the various values considered of u p , u z (i.e., the parameters setting the strength of the interface-adsorption forces orienting and keeping the NCs at the interface plane, respectively). The black line borders indicate the periodic boundary conditions. The graphs above each snapshot show, for each simulation and with respect to t=μs, the Lennard-Jones potential (light blue line) U LJ =10 −17 J summed over all the bead-bead pairs, the interface-adsorption potentials (green line) U z =10 −17 J and (fuchsia line) ðU φ þ U ψ Þ=10 −17 J, summed over all NCs, and the total internal potential energy (black line) U pot =10 −17 J; see Eq. (7).
FIG. 14. For the assembled-phase diagram of Fig. 5, in the case of total bonding energy between two NCs attached by f100g facets given by E ¼ −1.3 × 10 −19 J, we show here a top view of the interface at the end of each simulation (t ¼ 2 μs) for the various values considered of u p , u z (i.e., the parameters setting the strength of the interface-adsorption forces orienting and keeping the NCs at the interface plane, respectively). The black line borders indicate the periodic boundary conditions. The graphs above each snapshot show, for each simulation and with respect to t=μs, the Lennard-Jones potential (light blue line) U LJ =10 −17 J summed over all the bead-bead pairs, the interface-adsorption potentials (green line) U z =10 −17 J and (fuchsia line) ðU φ þ U ψ Þ=10 −17 J, summed over all NCs, and the total internal potential energy (black line) U pot =10 −17 J; see Eq. (7). a f111g facet parallel to the interface are too weak to matter, compared to the other forces. Indeed, in all the simulations presented in this work, the total potential U c , i.e., summed over all NC-NC pairs, remains always essentially constant, compared to the other potentials; hence, we do not show it in the energy plots. The total soft repulsive potential U R , i.e., summed over all bead-bead pairs, is also not shown in the energy plots of all the simulations presented in this work, since it remains essentially constant. Note, however, that U R plays a fundamental role since it prevents NC-NC attachment by f111g, f110g facets. In Figs. 12-14, we present a larger version of the assembled-phase diagram presented in Fig. 5, showing for each simulation performed, i.e., for each cell of Fig. 5, a snapshot of the final NC configuration, i.e., at t ¼ 2 μs, from a top view of the interface. Above each snapshot, we show for each simulation a plot, with respect to t, of (light blue line) the total energy of NC-NC attachment by f100g facets, i.e., U LJ summed over all bead-bead pairs, (green line) the total potential keeping the NCs at the interface plane, i.e., U z summed over all NCs, (fuchsia line) the total potential orienting the NCs with a f111g facet at the interface plane, i.e., U φ þ U ψ summed over all NCs, and (black line) the total potential internal energy U pot of the system; see Eq. (7).