Jessica Colombel
email: jessica.colombel@inria.fr

David Daney
email: david.daney@inria.fr

Vincent Bonnet
email: bonnet.vincent@gmail.com

François Charpillet
email: francois.charpillet@inria.fr

Markerless 3D Human Pose Tracking in the Wild with fusion of Multiple Depth Cameras: Comparative Experimental Study with Kinect 2 and 3

Human-robot interaction requires a robust estimate of human motion in real-time. This work presents a fusion algorithm for joint center positions tracking from multiple depth cameras to improve human motion analysis accuracy. The main contribution is the proposed algorithm based on body tracking measurements fusion with an extended Kalman filter and anthropomorphic constraints, independent of sensors. As an illustration of the use of this algorithm, this paper presents the direct comparison of joint center positions estimated with a reference stereophotogrammetric system and the ones estimated with the new Kinect 3 (Azure Kinect) sensor and its older version the Kinect 2 (Kinect for Windows). The experiment was made in two parts, one for each model of Kinect, by comparing raw and merging body tracking data of two sided Kinect with the proposed algorithm. The proposed approach improves body tracker data for Kinect 3 which has not the same characteristics as Kinect 2. This study shows also the importance of defining good heuristics to merge data depending on how the body tracking works. Thus, with proper heuristics, the joint center position estimates are improved by at least 14.6 %. Finally, we propose an additional comparison between Kinect 2 and Kinect 3 exhibiting the pros and cons of the two sensors.

Introduction

Human movement observation is necessary in robotics as for other applications to analyze human behaviours or to control assistance systems. When the perceptive system is used in situ (e.g., at home or at work) the usability of the system is almost as important as the minimum metrological performance to ensure the validity of the results obtained.

Nowadays stereophotogrammetric systems such as Vicon, Qualysis or Optitrack systems are used as golden standard for position tracking for laboratory experiments. However, these systems are too expensive and difficult to install on site. Furthermore, these may be constraining for the observed human because of the markers' placement. These constraints prevent the natural observation of humans in their environments such as home or work and do not permit to handle real life situations.

The emergence of wearable sensors or depth cameras capable of measuring human movement has extended the scope of investigation for personal monitoring outside the laboratory. Numerous algorithms based on Machine Learning for vision have obtained appealing results on human motion observation [START_REF]Kinect -Windows app development[END_REF][START_REF] Cao | OpenPose: Realtime Multi-Person 2D Pose Estimation using Part Affinity Fields[END_REF]. Nevertheless, the accuracy obtained and the reliability of skeleton tracker results are not always satisfactory for biomechanical analysis or to ensure safe interactions with assistance robots. However, the use of several devices and fusion techniques can deeply improve such systems with a good compromise between quality, usability and cost. In particular, it is interesting to favour data redundancy for visual systems that suffer from occlusions [START_REF] Moon | Multiple Kinect Sensor Fusion for Human Skeleton Tracking Using Kalman Filtering[END_REF].

In this context a Kalman filter based on a human model is proposed to merge and filter skeleton data from multiple sensors in real time. The Kalman filter is well known for the fusion of data from sensors, even from different types of sensors (i.e. depth camera and wearable sensors as IMU) [START_REF] Feng | Fusing Kinect sensor and inertial sensors with multi-rate Kalman filter[END_REF]. Especially for data from depth sensors, this type of algorithm can be also used with anthropometric constraints and a dynamical model for human motion to compensate for the lack of a model in the body tracking software [START_REF] Colombel | Physically Consistent Whole-Body Kinematics Assessment Based on an RGB-D Sensor. Application to Simple Rehabilitation Exercises[END_REF]. One objective of this paper is to evaluate if this approach remains useful and meaningful for the body tracking software Kinect 3 (Azure Kinect). Technical improvements over the body tracking software Kinect 2 are expected and are promising given the improvements observed between the first two generations of Kinect and the technical characteristics of both Kinect (See Table 1) [START_REF] Wang | Evaluation of Pose Tracking Accuracy in the First and Second Generations of Microsoft Kinect[END_REF][START_REF] Pagliari | Calibration of Kinect for Xbox One and Comparison between the Two Generations of Microsoft Sensors[END_REF].

In this paper we are comparing results of the proposed fusion algorithm on one experiment in two parts: one with Kinect 2 and the other with Kinect 3 body trackers (see Figure 1), with the Qualisys system as reference for the ground truth. The main contributions of this paper are:

• a skeleton fusion algorithm based on an Extended Kalman Filter (EKF) with human model independent of sensors; • a comparison of skeleton tracking of Kinect 2 and Kinect 3 faced to a reference system.

Human model and fusion methods

Joint center positions estimates

Both Kinect sensors provide skeleton estimate through the estimate of the 3D Joint Center Positions (JCPs). Kinect 2 estimates 25 JCPs including 3 joints for each hand, while Kinect 3 estimates the 26 JCPs including 5 joints for the head and 1 for each clavicle. Moreover, the joint center positions located on the spine are also different for each body tracker (namely the SpineChest and SpineNaval for Kinect 3 versus SpineShoulder and SpineMid for Kinect 2). Since both Kinects have different JCP estimates, we have selected 15 of them, called retained JCPs, located on the arms (except hands), legs and neck to be compared and analyzed. Figure 2a shows both skeleton models and the retained JCPs for the comparison. Other remaining joints were not considered in this analysis. The same 15 JCPs were also estimated using a reference stereophotogrammetric system and a set of 39 retroreflective-markers. These markers were set to match the Plug-In Gait full body template popularized by Vicon for the placement of retrofelective-markers [START_REF] Davis | A gait analysis data collection and reduction technique[END_REF]. Additional markers located on the lateral side of the knees, ankles and elbows were integrated to this template (see Figure 2b). Doing so reference JCPs were calculated using regression equations. For the elbows, wrists, knees and ankles it consisted to calculate the average position between medial and lateral markers. The hip JCPs were estimated using the Bell's regression equation [START_REF] Bell | Prediction of hip joint centre location from external landmarks[END_REF] based on the four pelvis markers. The shoulders JCPs were estimated by assuming an offset under the acromion markers along the vertical axis defined by the trunk segment [START_REF] Rab | A method for determination of upper extremity kinematics[END_REF]. Fig. 2b shows the location of the considered retroreflective markers and the estimated reference JCPs.

Extended Kalman Filter

The proposed EKF is based on a biomechanical model representing the human locomotor system. It is used to relate the JCPs estimated by the Kinect with the joint kinematics trough a kinematics model having fixed segment lengths. This kinematics model, composed of = 31 revolute joints ( ) and = 12 segments, was defined according to the recommendations of the International Society of Biomechanics [START_REF] Wu | ISB recommendation on definitions of joint coordinate system of various joints for the reporting of human joint motion-part I: ankle, hip, and spine. International Society of Biomechanics[END_REF][START_REF] Wu | ISB recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion-Part II: shoulder, elbow, wrist and hand[END_REF]. Modified Denavit-Hartenberg convention [START_REF] Khalil | SYMORO+: A system for the symbolic modelling of robots[END_REF] was used to calculate the forward kinematics model that allows to express the position of each of the = 18 estimated JCPs as a function of the joint angles ∈ R and of segment lengths L ∈ R . This forward kinematics model is used as the measurement model ℎ in the proposed Kalman filter. The EKF estimates the state vector X = [ L] while tracking the measurement vector Z ∈ R 3 composed of 3D coordinates of JCPs provided by each Kinect's skeleton data. The state and measurement vectors are modelled as follows at time :

X = (X -1 ) + w -1 , Z = ℎ(X ) + v , (1)
where is the state model, w represents system noise defined by (w) ∼ N (0, Q) with Q the model covariance noise matrix, and v represents the measurements noise defined by (v) ∼ N (0, R) with R the measurement covariance noise matrix.

The proposed state model is approximated by a linear form and denoted F. It assumes that between two consecutive samples joint angles and segment lengths are constants, which is suitable for slow motion. Moreover, the Jacobian matrix H defined by ℎ X is used as the local linearization of the measurement model. Thus, the equations of the extended Kalman filter are written as follows, with, at each time the prediction phase first:

X = F X -1 (2)
P = F P -1 F + Q (3)
and then the update phase :

K = P H (H P H + R) -1 (4) X = X + K (Z -ℎ(X )) (5) P = (I -K H ) P (6)
with P the error co-variance matrix, K the Kalman's gain matrix and I the identity matrix with the good size.

Sensors data fusion

There are two methods for merging multiple sensors with Kalman as a function of the type of fused variables [START_REF] Roecker | Comparison of two-sensor tracking methods based on state vector fusion and measurement fusion[END_REF]. They are called state-vector fusion for the merge of the state vector X and measurement fusion for the measurement vector Z. According to Gan et al [START_REF] Gan | Comparison of two measurement fusion methods for Kalman-filterbased multisensor data fusion[END_REF], these two methods are functionally equivalent when the ℎ measurement models of sensor are identical for each sensors to merge. However state-vector fusion requires to augment the measurement vector while measurement fusion merges observation before filtering data with the Kalman filter. Thus for state-vector fusion the computational cost grows with the number of sensors. As the proposed algorithm tends to operate in real time with a variable number of sensors and the measurement model ℎ are identical for all sensors, the measurement fusion was chosen. Thus the proposed method weights the measurements at each time without increasing the matrix sizes as follows:

Z =       =1 R -1       -1 =1 R -1 Z (7) H =       =1 R -1       -1 =1 R -1 H (8) R =       =1 R -1       -1 (9)
with Z the measurement vector, R is the covariance matrix of the measurement noise and H the Jacobian of the measurement model of sensor at time , and the number of sensors, here the number of Kinects.

Heuristics

The fusion of multiple sensors can be done without the knowledge of sensor confidence index. However, in the literature, several heuristics were used to weigh the measurements for merging multiple Kinect 1 or Kinect 2, such as point continuity with respect to velocity [START_REF] Yeung | Improved Skeleton Tracking by Duplex Kinects: A Practical Approach for Real-Time Applications[END_REF], segment length stability [START_REF] Moon | Multiple Kinect Sensor Fusion for Human Skeleton Tracking Using Kalman Filtering[END_REF] or placement with respect to Kinect [START_REF] Otto | Presenting a Holistic Framework for Scalable, Marker-less Motion Capturing: Skeletal Tracking Performance Analysis, Sensor Fusion Algorithms and Usage in Automotive Industry[END_REF]. These heuristics are based on technical specifications of the manufacturer depending on the sensors. So as to improve the fusion result, we propose a new set of heuristics inspired by those cited above that represents a confidence index in each sensor. Even though we merge several Kinect 2 on one side and several Kinect 3 on the other, they are based on the characteristics of the Kinect 2 only. These heuristics result in weights ∈]0.1]. Two of them are related to the placement of the body in relation to the sensor:

=              0.1 if < 1.5 1 if 1.5 ≤ < 3.5 0.5 if 3.5 ≤ < 4.5 0.1 else. (10) = -2 • + 1 (11)
with the depth distance of the spine base from the sensor and with the rotation of the body from the Kinect. The two following weights represent the confidence given in the sensor for the global body tracking trust (number of joints tracked over the total number of joints observed by the body tracking) and the time synchronisation to merge only closer frames (observations from different sensors must be close in time to be synchronized):

= (12) = -0.75 Δ • Δ + 1 ( 13
)
with and respectively the number of joints seen and the total number of joints of the camera as well as Δ and Δ the maximum time between two accepted frames and the time between two frames, respectively. The number of joints seen are known thanks to the tracking state variable given by the body tracker to know if the joint is measured or estimated. These values are then used to weight the R matrix of each sensor .

Calibration

Calibration of visual sensors is specially of importance when merging data coming from multiple sensors. In this study, solely the skeleton data coming from several Kinect sensors are merged. Skeleton data consist in a labelled point cloud representing the whole body at each time. The objective of the calibration process is to align these point clouds. Haralick et al. [START_REF] Haralick | Pose estimation from corresponding point data[END_REF] proposed a method to estimate a transformation matrix from 2 paired sets of = 1, . . . , 3D points:

Z 1 = [z 1,1 , . . . , z 1, ] and Z 2 = [z 2,1 , . . . , z 2, ] , with dim(z , ) = 3 × 1, = 1, 2. Let D = 1 =1 =1 (z 2, -z2 ) (z 1, -z1 ) (14) with z1 = =1 z 1, =1 and z2 = =1 z 2,

=1

Then {U, S, V} = SVD(D) [START_REF] Gan | Comparison of two measurement fusion methods for Kalman-filterbased multisensor data fusion[END_REF] such that D = U.S.V C = U.

1 0 0 0 1 0 0 0 det( . ) .V (16)
t = z1 -Cz 2 ( 17
)
where Z =1,2 are the measurement vectors composed of 3D JCPs for each sensors = 1, 2 and are weights of point (0 or 1 based on the tracking state variable given by the body tracking SDK). The rotation matrix C and the translation vector t represent the positioning of the sensors = 2 with respect to = 1.

A singular value decomposition (Eq. 15) is used to solve this orthogonal Procrustes problem (Eq. 16 ). To identify the 6 parameters of the transformation matrix in a robust way, the information provided by the point clouds is cumulated on several frames. The identification continues until the updating of the parameters converges to constant values. Unlike Kinect 2 which has a tracking state variable to give a quality index on each joint, Kinect 3 does not have this type of confidence index. Thus, we propose to use a median type method to select points for correspondence. For each frame, distance between each joint to the spine base are calculated for both sensor and compared. The 50% of points with a fewer error distance between sensor are taken ( = 1).

Once all the transformation matrices between cameras are calculated, the fused body can be express in a common base. This method is useful when it is needed to estimate the absolute position of the human body in the environment.

It should be noted that, to facilitate the correspondence between the three systems (multiple Kinect 2 and Kinect 3 and reference systems) and avoid error due to calibration, the common base chosen for this work is the skeleton itself. The base is defined by the spine base as the origin, the hip right-hip left axis for the x-axis, the spine axis for the y-axis and the z-axis orthogonal to the other two axes. To insure an orthonormal base, the Gram-Schmidt method was used.

Experiments and results

Experimental setup and protocol

Four tasks requiring whole body motions were considered in this study: squat with lateral extensions of arms, stepping in position with both legs, body tilt, and waving arms. Each task was performed with ten consecutive repetitions for each movement. All trials started and finished in a resting position: standing, arms along the body. These tasks were chosen to cover most degrees of freedom for arms, legs and trunk.

Whole-body motion was simultaneously collected by a reference stereophotogrammetric system consisting of 8 infrared cameras (Qualisys, 39 markers, 30 Hz) and by three Microsoft Kinect 2 (SDK, 30 Hz) and three Kinect Azure (NFOV, SDK v. 9.2, 30 Hz). The Kinect sensors were located in a triangle fashion with two Kinects in front (0°) of the participant, two on his right side (80°) and two on his left (-80°) side. Fig. 3 shows an overview of the experimental setup. Each Kinect was connected to a Windows 10 computer and all data were streamed on a ROS-based software to be elaborated on-line. Kinect 3 were synchronized with 3.5 mm audio cables in a daisy chain configuration. Due to interference between Kinect 3 and the reference system, the master cable has been split to trigger the reference system with the signal of the Kinect 3 master. Moreover, Kinect 3 body tracking SDK is badly disturbed by the motion capture retro-reflective markers. To minimise theses disturbances miniature markers of 2.5 mm of diameter were used.

Performance analysis

Kinect 2 has shown satisfactory results when the human body is facing the sensor without occlusion. Thus, to highlight the contribution of the proposed fusion algorithm, only the side sensors were considered (-80°and 80°), with heuristics and without heuristics, now referred to as Fusion H and Fusion !H, respectively. This was done for both versions of the Kinect sensor. The ability of the proposed fusion algorithm was then assessed by calculating the accuracy through the Euclidean distance between filtered JCPs and those obtained from the reference stereophotogrammetric system (see section II.A). The accuracy of the raw data gathered from the three sensor positions was also investigated for both Kinects. This was done to better understand the limitations of their specific body tracking software. Moreover, the mean and standard deviation of segment lengths were calculated on the raw and filtered data to attest to the consideration of anthropomorphic constraints by reducing the variability on segment lengths.

Results

Joint center position estimate accuracy assessment

First, the accuracy of the JCP estimates is evaluated by calculating the Euclidean distance between the JCPs obtained from the reference system and the JCPs estimated with the Kinect sensors. This accuracy is an absolute error in position. Table 2 reports mean and standard deviation distance for each joint. This is calculated for both types of Kinects when using raw data with different sensor's positions (0°, -80°and 80°) and when using merged data with and without heuristics.

As expected for Kinect 2, the Fusion H shows a better accuracy with a slightly shorter average distance (106 ± 41 mm) than Fusion !H (114 ± 36 mm). Fusion H improves on average of 14.6% for left (-80) and 28.9% for right (80) Kinect 2 results. When Fusion !H is used the improvement is of 8.0% and 23.4% for the left and right sensor, respectively. Globally, fusion algorithms tend to improve each joint accuracy and seem to be slightly equivalent. However the left elbow and wrist present a large difference between Fusion H and Fusion !H with both 121 mm for Fusion H, and 201 and 183 mm for Fusion !H, respectively. For these two joints, the results are particularly different between the right and left sensors, more than twice the mean and standard deviation. The difference on these two joints leads to the global difference between placements. It can be explain by the fact that the motion was not exactly symmetrical and the self-occlusion affected more these joints for the right sensor than the other.

Unlike its previous version, Kinect 3 shows better results when using Fusion !H (84 ± 36 mm) than when using Fusion H (99 ± 206 mm). Both side Kinect 3 results are improved by Fusion !H with 21.4% and 16.8% for the left and right Kinect 3, respectively. The improvement of raw data by Fusion H are less good with 7.6% for the left sensor and only 2.2% for the right sensor. Finally, the accuracy of the feet position obtained when Fusion H algorithm is used decreases greatly. This can be observed with the very high standard deviations of the feet (1187 mm and 1259 mm).

Segment length variation

There is a consensus in the community stating that the body tracking of Kinect 2 is not anthropometry-based, resulting in segment length variations. Thus, numerous researchers tried to prevent segment length variations to improve human body tracking [START_REF] Tripathy | Constrained Particle Filter for Improving Kinect Based Measurements[END_REF][START_REF] Skals | A musculoskeletal model driven by dual Microsoft Kinect Sensor data[END_REF]. The proposed Extended Kalman Filter used the process matrix covariance Q to force the length converge to a constant value. Fig. 4 represents the dispersion over all trials of the segment length. It shows that segment lengths of Kinect 2 have greater dispersion and more outliers than Kinect 3 for raw data. Moreover it illustrates the constraints on segment length achieved with the proposed Kalman filter. The segment length stability improves significantly with the Kalman filter. This is visible on the comparison between raw and filtered data for the right lower leg. There is also a strong reduction of the outliers. Table 3 confirmed these observations with an average standard deviation obtained for Kinect 2 of 19 mm instead of 8 mm or 9 mm when filtering.

It is interesting to notice that the Kalman filter also slightly improves the accuracy of the average absolute error of segment length estimate (37 mm for raw data instead of 35 mm for Fusion H). On the contrary, results obtained with the Kinect 3 show that Fusion !H degrades length estimate with an absolute error of 20 mm instead of 17 mm without filtering. Finally, Fusion H has terrible results on standard deviation, particularly for feet.

Additional results

These experimentations with both sensors raised some additional results that highlighted the different characteristics between Kinect 3 and Kinect 2. As mentionned in section 3.2 and confirmed by the results in Table 2, the central Kinect 2 shows much higher accuracy than those located on the side (82 ± 36 mm versus 124 ± 54 and 149 ± 71 mm). It is interesting to note that the Kinect 3 presents opposite results. The accuracy is lower for the central Kinect (121 ± 85 mm) than for the ones located on the side (107 ± 59 and 101 ± 53 mm).

Discussion and limits

The proposed fusion algorithms show improved accuracy on both Kinect 2 and Kinect 3. The differences due to the choice of heuristics were expected since the two Kinects do not present the same technical characteristics, e.g. the field of view and the depth resolution. Moreover, Kinect 2 was designed to view people facing the sensor for video games, but Kinect 3 is presented as an industrial sensor with a body tracking algorithm designed for multiple purposes, including back and front recognition. Several other differences as directly-included time synchronization for Kinect 3 and the absence of a tracking state variable for each joint makes at least two of the heuristics unnecessary. It can be concluded that the proposed algorithm improves body tracking but can still be enhanced by using sensor specific heuristics.

As expected, thanks to the use of the Kalman filter and anthropomorphic constraints, the length estimates vary less, although their estimation may be less accurate. The proposed approach is really interesting for Kinect 2 but seems to be less useful for Kinect 3. Unlike Kinect 2, Kinect 3 promises to provide an anatomically consistent skeleton tracking.

Despite these improvements on Kinect 3, it does not appear to greatly improve the accuracy of skeletal observation. The results obtained with the central Kinect 3 are quite surprising given that they are particularly poor compared to the ones obtained with the side sensors. However, even if the body tracker improves, the fusion will still be necessary to counter the effect of occlusion. These experiments gave additional results on the new Kinect 3 such as interference with stereophotogrammetric systems and reflective markers. The use of miniature markers of 2.5 mm of diameter reduces the effect of interference with the Kinect 3. In addition, the synchronization between the sensors constrains the Qualisys system to work at 30 Hz instead of 300 Hz. Naembadi et al. studied this question with Kinect 2, but Kinect 3 seems to be much more sensitive to interference [START_REF] Naeemabadi | Investigating the impact of a motion capture system on Microsoft Kinect v2 recordings: A caution for using the technologies together[END_REF]. Another point to highlight is the computational cost for running Kinect 3. It requires a good GPU (NVIDIA GEFORCE GTX 1070 or better) to work at 30 Hz. These results could question the portability of the sensor for embedded systems. It is interesting to notice that, if the portability of the system is critical and that the visual system as depth camera is not adapted to the application (i.e. the person had to be observed when moving in large area), other approach can be used [START_REF] Banos | Kinect=IMU? Learning MIMO Signal Mappings to Automatically Translate Activity Recognition Systems across Sensor Modalities[END_REF].

Conclusion

This paper has presented a fusion algorithm based on an Extended Kalman Filter for merging skeleton data obtained from multiple depth cameras. The fusion algorithm was tested with the new Kinect 3 and with the Kinect 2 with regard to a reference provided by a stereophotogrammetric system. The results showed that this algorithm is robust to a change of body tracker software and overall improves joint center positions estimate. When the proposed fusion algorithm is used, the JCP estimate accuracy improved by at least 14.6 %. The results prove that with or without fusion heuristics, the algorithm improves accuracy for each type of sensor. As shown for Kinect 2, a good definition of the heuristics allows even better results. However, the heuristics chosen to improve the fusion algorithm of multiple Kinect 2 are not suitable when using Kinect 3 sensors. This is due to differences between both sensors such as time synchronization, the lack of tracking confidence index number in Kinect 3 or to the fact that Kinect 3 is able to recognize back-front positioning. In addition to this last point, which is very important for usability, the improvement of segment length estimates is the main improvement obtained when the proposed algorithm is applied to Kinect 3 rather than Kinect 2. These two improvements open the door to a better observation of humans in their environment with the possibility of observing them at 360°.

Fig. 1

 1 Fig. 1 Pictures of the last two generations of Kinect : (Top) Kinect v2, (Bottom) Kinect Azure

Fig. 2 :

 2 Fig. 2: Location of (a) joint center positions of the Kinect 2 (yellow), and of the Kinect 3 (light blue). The green circles represent the 15 retained JCPs that stand for common JCPs to compare. (b) Retroreflective markers used for the reference stereophotogrammetric system (dark blue) and the estimated reference joint center position (green)

Fig. 3 :

 3 Fig. 3: (a) Schematic top view of the experimental setup: black circle is the human placement, black camera represents a Kinect 2 and blue camera a Kinect 3. (b) Panorama showing the location of the Kinect sensors facing the subject.

Fig. 4 :

 4 Fig. 4: Boxplot representing the dispersion of length on all trials for Kinect 2 and Kinect 3 with raw data and filtered data. Raw Kinect 2 is dark red, Fusion H Kinect 2 is orange, Fusion !H Kinect 2 is yellow, raw Kinect 3 is blue, Fusion H Kinect 3 is light blue and Fusion !H Kinect 3 teal.

Table 1 :

 1 Table of technical characteristics of Kinect 2 and 3.

		Kinect 2		Kinect 3 Azure
			NFOV NFOV binned WFOV WFOV binned
	Depth Camera (pixel) 512×424 640x576 320x288 1024x1024	512x512
	FPS	30	30		15	30
	Min depth distance (m)	0.5	0.5		0.25
	Max depth distance (m)	4.5	3.86	5.46	2.21	2.88
	Horizontal FOV (degree)	70	75			120
	Vertical Fov (degree)	60	65			120
	Skeleton joint define	25		26 (SDK v.9.2)

Table 2 :

 2 Joint center positions estimate obtained with the Kinect 2 and the Kinect 3 for raw (0°, -80°and 80°) and filtered data from the fusion algorithm with (Fusion H) and without heuristics (Fusion !H). Accuracy as absolute error is reported as Mean ± Standard Deviation [mm] over all the analyzed trials. Results in bold present the best accuracy for each joint between Fusion algorithms and raw data (except the facing camera (0°)).

	Kinect 2
	Joint name

Table 3 :

 3 Segment length estimation obtained by Kinect 2 and Kinect 3 for raw data (independently of the placement) and filtered data from fusion algorithm with (Fusion H) and without heuristics (Fusion !H) comparing to segment length reference. Accuracy as absolute error of segment length of interest has been reported as Mean ± Standard Deviation [mm] over all the analyzed trials. Upper Leg R 85 ± 35 83 ± 11 87 ± 12 42 ± 15 66 ± 67 57 ± 5 Lower Leg R 41 ± 34 30 ± 14 31 ± 11 22 ± 12 31 ± 84 9 ± 3 Foot R 53 ± 13 59 ± 14 62 ± 10 5 ± 4 100 ± 4275 9 ± 4 Upper Leg L 94 ± 31 82 ± 11 86 ± 12 40 ± 15 66 ± 67 56 ± 5 Lower Leg L 24 ± 22 29 ± 14 30 ± 10 20 ± 12

		Joint name			Kinect 2				Kinect 3
						raw	Fusion H Fusion !H	raw	Fusion H Fusion !H
		Clavicle R	13 ± 15	9 ± 4		8 ± 4	5 ± 4	10 ± 36	3 ± 2
		Upper Arm R 33 ± 18 53 ± 8		58 ± 9	9 ± 7	19 ± 25	19 ± 5
		Arm R	14 ± 12 7 ± 6		9 ± 7		8 ± 6	20 ± 24	23 ± 8
		Clavicle L	16 ± 19	9 ± 4		8 ± 4	8 ± 5	10 ± 35	3 ± 2
		Upper Arm L 31 ± 16 52 ± 8		57 ± 9 10 ± 8 18 ± 25	18 ± 5
		Arm L	14 ± 14 7 ± 6		9 ± 7		9 ± 7	20 ± 24	23 ± 8
		Hip R	18 ± 10	4 ± 2		3 ± 2 21 ± 3 40 ± 73	23 ± 1
		Hip L	25 ± 10	4 ± 2		3 ± 2	30 ± 3	40 ± 73 23 ± 1
		0.4	K2 raw K2 Fusion H K2 Fusion !H K3 raw K3 Fusion H K3 Fusion !H
	Length (m)	0.2
		0
		Clavicle	Upper	Arm	Clavicle	Upper	Arm	Hip	Upper	Lower	Foot R	Hip	Upper	Lower	Foot L
		R		arm R	R	L	arm L	L	R	leg R	leg R		L	leg L	leg L
									Segment

32 ± 84 10 ± 3 Foot L 53 ± 14 57 ± 9 62 ± 9 15 ± 7 108 ± 4650 7 ± 4 Total 37 ± 19 35 ± 8 37 ± 9 17 ± 8 41 ± 691 20 ± 4

Jessica Colombel, David Daney, Vincent Bonnet, and François Charpillet

Acknowledgement

This work has been partly funded by the CPER IT2MP (Contrat Plan État Région, Innovations Technologiques, Modélisation Médecine Personnalisée) and FEDER (Fonds européen de développement régional).


Michael Thomas

Irene Suarez-Martinez

Li-Juan Yu

Amir Karton

Graham Chandler

Marc Robinson

Isabelle Cherchneff

Dahbia Talbi
email: dahbia.talbi@umontpellier.fr

Dino Spagnoli
email: dino.spagnoli@uwa.edu.au

Irene Suarez Martinez

Atomistic simulations of the aggregation of small aromatic molecules in homogenous and heterogenous mixtures

come   L'archive ouverte pluridisciplinaire

Introduction

Weak intermolecular forces control the physical properties of aromatic molecules and play a crucial role in processes such as protein-ligand recognition, 1 the packing and assembly of aromatic molecules in crystals, 2 and soot formation. [3][4][5] These interactions between aromatic molecules are strong enough so clusters can form via intermolecular forces. 5 The early chemistry of soot formation is now generally thought to be controlled by these nonbonding interactions. Understanding the formation of soot by aromatic molecules has received a great deal of research mainly due to the environmental implications. 3 Aggregation of aromatic molecules is also of special interest in the field of astrochemistry, in particular for the study of the chemical composition of the outer sphere (photosphere) of carbon-rich stars. Evidence has been accumulating for years on the similarity between the chemical pathways observed in combustion processes on Earth, which form benzene, polycyclic aromatic hydrocarbons (PAHs) and soot, with the high temperature conditions in circumstellar regions. [6][7][8] Observation from radio telescopes have determined that an envelope of carbon dust surround red-giants. These particular stars emit more carbon and less oxygen than most other stars and are the main carbon dust providers to local galaxies. 9 Theoretical studies of the growth of PAHs 6,7,10,11 is commonly explained through hydrogen abstraction and acetylene addition (i.e., HACA mechanism 7 ). PAHs could form and grow in a distinctive temperature range of 900-1500K, but the yields were rather low to explain the ubiquity of dust in carbon stars. The mechanism of aggregation of PAHs has not yet been addressed, and it is key to understanding the production of carbon grains from stars.

Similarities exist between the chemical and physical conditions of acetylene-rich flames and those pertaining to the upper atmosphere of carbon stars where dust forms. Both environments are rich in hydrogen, carbon, acetylene and oxygen, and have similar temperature regimes. Experiments in the laboratory on carbon dust formation mimicking the low-density, warm environments of carbon star winds have been undertaken. 12,13 Plasma discharge experiments have shown that the growth of unsaturated hydrocarbons and aromatics proceeded through chain growth and ring formation. 14 More recently, the production of carbon clusters in ultra-high-vacuum gas aggregation sources shows that carbon growth mainly proceeds through the condensation of pure carbon chains and rings, where aromatics represent a minor population of carbon clusters. 15 Theoretical studies related to soot formation shows that the aggregation of PAHs through London dispersion interactions can be efficient for a large range of temperature from 400 K to 1600K. 16 To advance our understanding of PAH aggregation into small carbon clusters, we conducted this extensive and systematic investigation through molecular dynamics (MD) and ab-initio simulations. We considered both homogenous (one molecular species) and heterogeneous (multi-molecular species) aggregation between 100K and 500K of three small aromatic molecules: benzene, naphthalene and anthracene. The literature on aggregation of polyaromatic hydrocarbons via simulations is extensive. [16][17][18][19] Moreover, the mechanism and the molecules involved in soot nucleation is unknown. 4,20 This study provides insight into the mechanism of aggregation of small aromatic molecules. We believe that this is the first time that a systematic investigation of heterogenous mixtures of small aromatic molecules has been performed. In-depth analysis arrangements of the molecules within the cluster is performed and compared to ab-initio calculations of dimers and trimers. Finally, our results are discussed within the astrophysical context of dust formation in evolved stars.

Methods

Molecular Dynamics Simulation Details

Seven systems were constructed for the MD simulations. The composition of each system is detailed in Table 1. The number of molecules in each system was chosen to maintain a constant atomic density using periodic boundary conditions in a cubic cell of length 120 Å. The size of the simulation cell and number of molecules were chosen because we wanted to simulate a system where aggregation would occur on a timescale accessible to MD. We note that the density of the simulations is far higher than densities found in space. All systems were simulated with the molecular dynamics package NAMD 2.12 21 . All simulations were conducted using the NVT ensemble (constant number of atoms, temperature and volume) with the Langevin thermostat 22,23 and a dampening parameter of 5 ps to maintain constant temperature. Particle Mesh Ewald 24 was used to calculate electrostatic interactions. A timestep of 1 fs was used in all cases.

Initially the OPLS-AA force field 25 was chosen based on previously reports by Takeuchi who calculated the potential energy of PAHs clusters in static simulations using this force field 26,27 . However, our dynamic simulations showed that aggregation was only observed when benzene was present using the OPLS-AA force field. Simulations involving naphthalene and anthracene did not show aggregation even at temperatures below the condensation point. We repeated the calculations using the parameters from CGenFF within the CHARMM36 force field 28 to describe bonding and non-bonding interactions and aggregation occurs for all systems. All results that follows use the CHARMM force field.

All systems were modelled for five temperatures: 100, 200, 300, 400 and 500 K. Each temperature was modelled five times using five different initial configurations, which is a method used in other MD studies. 29 Initial configurations were produced from randomly positioned molecules heated at 1000 K for 100, 200, 300, 400 and 500 ps. This method ensures the final configurations are not biased by the initial configuration. All simulations were performed for 10 ns. All seven systems displayed aggregation at 200 K, therefore these simulations were run until a single cluster was formed at this temperature, typically 100-200 ns. Each simulation trajectory was analysed to determine cluster sizes by searching atoms from different molecules within 4 Å. This distance was chosen as literature has shown that the maximum distance between two aromatic rings in a parallel displaced orientation would be no more than 3.9 Å. 30,31

Ab initio Simulation Details

High-level double-hybrid density functional theory (DHDFT) is used to investigate the geometries of dimer and trimer complexes. All geometry optimization, frequency, and energy calculations were performed using the Gaussian 09 program suite. [START_REF] Frisch | Gaussian 09[END_REF] In the first approach all the degrees of freedom have been fully optimized at the M06-L/6-31G(2df,p) level of theory [START_REF] Zhao | [END_REF]34 in order to obtain the local minima structures on the potential energy surface (PES). We were able to apply this approach for the benzenebenzene, benzene-naphthalene, and benzene-anthracene dimers in the parallel displaced (PD) configuration. In all cases harmonic vibrational frequencies have been calculated to confirm each stationary point as an equilibrium structure (i.e., all real frequencies). In the second approach, we use the optimized monomers and scan the vertical intermolecular distance. The DHDFT PWPB95-D3 functional in conjunction with the Def2-QZVPP basis set 35,36 was used to carry out the PES scans. Empirical D3 dispersion corrections 37,38 are included using the Becke-Johnson 39 damping potential as recommended by Grimme et. al. 40 (denoted by the suffix D3). The Resolution of Identity (RI) approximation was used to speed up the calculations. This level of theory has been shown to give very accurate interaction energies for weakly bonded systems. 36 This approach is used to obtain the structures of all the dimers and trimers considered in the present work. It is worthwhile comparing the optimal intramolecular distances obtained via scanning only one coordinate (i.e., the R and R1 coordinates for T-and PD-structures, respectively) with those obtained via the fully optimizations. These results are gathered in Table S3 of the Supporting Information. The ring-to-ring distance R1 for fully optimized PD benzene dimer is 3.450 Å, whereas scanning this distance results in R1 = 3.620 Å. The difference between these two values amounts to 0.17 Å, that is 4.9% of the intramolecular bond distance. Similarly, the partial optimization leads to a bond length that is longer by 0.19 Å for the benzene-naphthalene dimer (i.e., 5.4% of the bond distance). For the benzene-anthracene dimer, the partial and full optimizations result in a similar intramolecular bond separation (the difference between the bonds being only 0.04 Å).

Results and Discussion

Clustering as function of temperature and mixture composition

The seven molecular systems were simulated at five different temperatures for 10 ns. The purpose of these simulations was to understand at which temperatures aggregation occurs within a practical MD timeframe for the different molecular systems. The maximum cluster size as a function of size is shown in the left-hand side of figures 1a-g. At 100 and 200 K, all molecular systems show evidences of cluster formation. For simulations involving only benzene there is very little aggregation at temperatures at 300 K and above. In the supplementary material table S1 shows all average number of clusters, average maximum cluster size and average cluster size for all replicates at each temperature. For simulations involving naphthalene and anthracene, clusters begin to form and are stable at higher temperatures of 300 K and 400 K. There is very little evidence of any clusters of more than two molecules forming for any molecular system at 500 K. For the same temperature, we find larger clusters in the heterogeneous mixtures compared to homogenous mixtures. For example, in the molecular system involving only benzene there is very little evidence of clusters forming at 300 K. However, in the heterogeneous molecular systems involving benzene and/or the polyaromatic molecules there is evidence of clusters forming at 300 K. Since the dominant interaction between these non-polar molecules is London dispersion forces, with its strength proportional to molecular weight, it is unsurprising that this would be the case. At 10 ns, the benzene homogenous mixture shows nearly half of the initial molecules aggregated into a single cluster at 200 K. There is very little evidence of any clusters forming above the 2-molecule threshold at 300, 400 or 500 K for the pure benzene molecular system. For the homogenous system of naphthalene, the bigger clusters are formed with increasing temperature from 100 K to 300 K, and no evidence of any cluster forming above 400 K (figure 1b). In the case of anthracene, clusters are observed at all temperatures except at 500 K (figure 1c). The calculation of the enthalpy change of sublimation from experimental studies for benzene ranges between temperatures of 193 41 -279 K, 42 naphthalene from 274 -353 K, 43 and anthracene from 337 44 -488 K. 45 Our results are in agreement with those values, and show that at temperatures above the sublimation point there is no aggregation. At temperatures below the point of sublimation, for each homogenous molecular system, aggregation occurs.

Intermolecular forces are stronger in heterogenous mixtures compared to homogeneous mixtures. Heterogenous clusters are observed at temperatures in which homogeneous mixtures showed no aggregation. For example, in the simulation involving benzene and naphthene mixture, clusters can be seen at temperatures up to 300 K (Figure 1d), whereas in pure benzene mixtures clusters only formed for temperatures under 200 K (Figure 1a). Similarly, clusters are formed in benzene and anthracene mixtures at temperatures up to 300 K (Figure 1e). In the mixture containing the two largest aromatic hydrocarbons (naphthalene and anthracene, Figure 1f), the biggest clusters are observed at 200 and 300 K, while at low temperatures (100 K) and high temperature (above 400 K), aggregation is none or limited. Interestingly, the mixture of the three aromatics do not form cluster above 400 K (Figure 1g), however pure mixtures of anthracene clearly formed clusters at that temperature.

The final configuration of the largest clusters formed at the end of 10 ns simulation is shown in the right-hand side of Figure 1 a-g. In heterogenous mixtures, clusters rarely form of one single molecular type. Moreover, larger aromatic molecules are often found in the centre of the clusters surrounded by smaller molecules. One particular exception is the mixture of the three aromatic molecules, in which at the 10 ns stage no pattern appears evident and cluster are often randomly arranged.

Molecular dynamics simulations at 200 K for > 100 ns

To investigate the final configuration and geometries of homogeneous and heterogeneous mixtures over longer time scale, the seven systems were simulated at 200 K. The temperature of 200 K was chosen because all molecular systems showed clustering at this temperature based on the 10 ns simulations. All replica simulations were run until all molecules in the calculations clustered into a single particle.

In the supplementary material table S2 shows the final MD time for the system to converge to form a single particle. Figure 2a shows the average cluster size across the five replicas for each mixture as a function of time. In the supplementary material figures S1-S7 show the variation of average cluster size, maximum cluster size and number of clusters for each of the replicates in each molecular system. All values in Figure 2a have been normalised by the total number of molecules in the system, therefore the value of 1 indicates a single cluster. first phase, a large number of small clusters quickly form within the first 1 ns. This is demonstrated by the small cluster size (Figure 2a) and large number of clusters (Figure 2b) at a time of 1 ns. In the second phase, after 1 ns, those smaller clusters aggregate, increasing the number of molecules comprising the cluster significantly (Figure 2a) and decreasing the number of clusters (Figure 2b). The alternative mechanism would display a constant small number of clusters, resulting in a much flatter curve for Figure 2b.

Structural analysis of clusters at 200 K

We analyse the arrangement of molecules within the clusters by pairs from the final configurations of the converged 200 K simulations. Three variables were used to describe the interactions between molecules: the distance between the centre of mass of the two molecules (r), and two angles (φ and θ) to describe the orientation with the two molecules with respect to each other. One angle is defined by the vectors normal to the ring plane of two molecules (φ), while the other (θ) angle is defined by two long axes of the molecules (see Figure 3a). The long axis in benzene is defined as the vector between the centre of mass and a carbon atom, while in naphthalene and anthracene it is defined through the longest length of the molecule. Due to the plane of symmetry through the ring plane, the range of φ and θ is 0 to 180°. The values plotted in the histogram have been corrected for sampling bias with a correction of 1/r 2 for r, and 1/sinφ and 1/sinθ for φ and θ respectively. The histograms have been coloured to illustrate angles between molecules at a large frequency (red), medium frequency (orange to yellow) and no frequency (black).

Figure 3b-g shows the structural analysis of molecule-molecule orientation for the entire length of the simulation. In the case of the benzene mixtures, molecules are between 4-6 Å apart and there is a large occurrence of configurations between the molecules in a T-configuration, φ=90°, and parallel-displaced (PD) configuration, φ=0° (see Figure 3b). Due to the symmetrical nature of benzene and rotation of the molecules, all possible θ values are observed with a slightly higher preference to 0°, 60°, 120° and 180° (Figure 3c). The most predominant configuration is more selective as the size of the PAHs increases.

In the case of the naphthalene-naphthalene interactions, the perfect T-configuration is observed less frequently than in benzene as depicted by the yellow region at φ=90°, instead we find a configuration localised around a tilted configuration at φ=45° and 135° (Figure 3d). The long axes of the naphthalene molecules are almost aligned parallel with each other (θ=0° and 180°) as seen in Figure 3. Anthraceneanthracene interactions display the only two distinct configurations with respect to the φ angle (Figure 3f): PD configuration (at φ=0°) and tilted configuration (at φ=45°). As shown by the black region at φ=90° in Figure 3f, the T-configuration is not found in any anthracene cluster. The arrangement of naphthalene and anthracene at φ=45° and θ=0° observed in the clusters is consistent with herringbone pattern found in the crystal structure of anthracene 46 (Figure 4). Heterogenous mixtures enables us to investigate the effect of another molecule in the structure of the cluster. In this paper we include r/φ analysis, while analysis of r/θ can be found in the supplementary material, figures S8-S9. Figure 5 shows the analysis for r/φ for all mixtures of two aromatic types.

Moving down the first column of Figure 5 we see the benzene-benzene interaction in the presence of benzene, naphthalene and anthracene. They are very similar to one another, indicating that these interactions are identical in the one and two component mixtures. The same is also true for the anthracene-anthracene interaction in the third column. The benzene-benzene and anthracene-anthracene interactions are not significantly affected by the presence of other species. However, the naphthalenenaphthalene interaction behaves differently depending on which other molecules are in the mixture. The most significant change in the naphthalene-naphthalene interactions occurs in the presence of anthracene. Anthracene is coercing the naphthalene-naphthalene interaction to behave more like the anthracene-anthracene interaction, thereby restricting the rotation, removing any T-shaped type configurations and producing more herringbone pattern within the cluster. Figure 6 shows the interactions in the three-component mixture. The benzene-benzene interaction and the anthracene-anthracene are the same as in the pure mixture and the two-component mixtures. Once again, the naphthalene-naphthalene interaction is modified by the presence of anthracene as in the anthracene and naphthalene mixture. Benzene-naphthalene interactions can rotate between T and PD configuration with no particular preference. Benzene-anthracene are dominated by tilted and PD configurations with limited T-configurations. Naphthalene-anthracene shows tilted and PD configurations with a total absence of the T-configuration, showing that anthracene restricts any Tarrangement.

Benzene-Benzene interaction

Anthracene-Anthracene interaction

Naphthalene-Naphthalene interaction

In the presence of: The final structure of each replicate of the three-component system from the longer simulations are displayed in Figure 7. All five final structures show a 'core' of anthracene in the centre of the cluster. Even when two smaller clusters merged, the anthracene cores in each cluster are seen to merge to form a new single core. As the interactions amongst aromatic molecules are primarily driven through the London dispersion forces, the strongest interactions will be the anthracene-anthracene interaction, and so the formation of the anthracene core aids to minimise the energy of each system. Naphthalene molecules fulfil the role of extending this anthracene core crystal. The naphthalene-anthracene interaction is weaker than the anthracene-anthracene interaction, and so naphthalene molecules are more mobile within the cluster than anthracene molecules. Benzene molecules provide the weakest interaction and remain mainly at the surface of the cluster.

where ∆𝐸 # is the binding energy of the dimer, 𝐸 # (𝐴) is the electronic energy of first molecule, 𝐸 # (𝐵) is the electronic energy of the second molecule and 𝐸 # (𝐴𝐵) is the electronic energy of the complex. A similar calculation is performed for the trimers, however, the number of molecules in the equation would be adjusted accordingly. The equilibrium intermolecular distances and binding energies are shown in Figure 8 as well as in Table S3 in the Supporting Information. The binding energy for dimers with benzene in PD configuration increases from 11.9 kJ mol -1 in the case of benzene-benzene to 19.3 and 23.9 kJ mol -1 in the case of the benzene-naphthalene and benzene-anthracene, respectively (Figure 8a). We note that our results for the benzene-benzene PD dimer are in excellent agreement with estimated CCSD(T)/CBS. In the case of trimer complexes in PD configuration, we obtain binding energies of 24.5 for the benzenebenzene-benzene, 39.7 for the benzene-naphthalene-benzene, and 48.9 kJ mol -1 for the benzeneanthracene-benzene. Thus, similarly to the PD dimers the binding energy for the benzene-anthracenebenzene trimer is twice as large as that in the benzene trimer, showing that benzene prefers aggregation with other aromatics rather than itself. The increase in the binding energy in both the dimers and the trimers is also associated with shorter distances (see Figure 8). It should be noted that when moving from the dimers to the trimers the binding energy per bond (or per p-p interaction) slightly increases. This trends also apply to the T-shaped dimers and trimers. Table S3 in the supplementary material shows all binding energies for the PD and T-shaped configurations at equilibrium distances for the dimers and trimers shown in Figure 8 as well as for dimers of naphthalene-naphthalene and naphthalene-anthracene.

Our high-level DHDFT results show a significant and systematic increase in the p-p binding energies when benzene is complexed with benzene, naphthalene, and anthracene. The same trends are observed for the PD dimers and trimers. These results are consistent with the MD results which show that anthracene tends to aggregate in the centre of the aromatic clusters. Our ab initio results show that Tconfigurations are less stable than PD. This is consistent with our MD simulations where Tconfigurations are less frequent than PD.

Conclusions

In the present study, we investigate the aggregation of different aromatic molecules at temperatures similar to those found in the dust formation zone around stars. Regions close to the star (1300 K) are rich in hydrogen and acetylene 50 promoting the formation of small PAHs. While the mechanism of formation of aromatic molecules is under debate, one possible mechanism is the growth of benzene to naphthalene, anthracene and larger PAHs can be pictured via a sequence of HACA processes. 7 Dust formation from carbon stars may thus result from a first phase of pure chemical growth of aromatics in the high-temperature, innermost regions followed by their aggregation and further growth once the gas temperatures have dropped to warm values (below 500K). Such a scenario leads to a population of small carbon clusters present at few stellar radii from the stellar surface.

The aggregation of PAHs into dust grains has not been considered in previous studies. The present study indicates the efficient clustering of small aromatic molecules through London dispersion forces at the warm temperatures (~ 300-500K) usually found at few stellar radii from the stellar surface. In our simulations, aggregation occurred for all systems containing larger polyaromatic molecules in this temperature range. Heterogeneous clustering will prevail at higher temperatures over homogeneous clustering, thereby facilitating the growth of non-planar aromatic clusters of small sizes, which may act as condensation nuclei. The structural analysis of our simulations reveals that benzene will orient itself in T-shaped and PD configurations. The structural behaviour of benzene molecules is due in part to the high degrees of symmetry within the molecule and is also further evidenced by the similar interaction energies between the two configurations. The polyaromatic molecules (naphthalene and anthracene) will not form a T-shaped configuration between molecules and will always form configurations localised around φ=45 and 135°, which is similar to the crystal packing of these systems as a solid. Moreover, in the molecular systems with all three molecules we observe evidence of a core of larger anthracene molecules surrounded by the smaller molecules. The intermolecular forces between anthracene molecules is the driving force to keep these molecules closer together in configurations close to the parallel displaced orientation.

The simulations presented in this study are the starting point for systematic and rigorous exploration of the aggregation of other aromatic molecules in homogenous and heterogenous mixtures. Recent laboratory studies, highlights the importance of small carbon chains in the dust formation process. 15 The role of carbon chains and the interaction with aromatic molecules in the production of small clusters in the dust formation zone of carbon stars needs further modelling investigation.

Figure 1 :

 1 Figure 1: Maximum cluster size averaged over the five replicas at 100 K (red), 200 K (green), 300 K (blue), 400 K (pink) and 500 K (cyan) for a) 60 benzene molecules, b) 40 naphthalene molecules, c) 30 anthracene molecules, d) 30 benzene and 20 naphthalene molecules, e) 30 benzene and 15 anthracene molecules, f) 20 naphthalene and 15 anthracene molecules, and g) 20 benzene, 13 naphthalene and 10 anthracene molecules. An example of the final configuration of the largest cluster formed for each molecular system is shown on the right of each plot.

Figure 2 :

 2 Figure 2: a) Average cluster size normalised by the total number of molecules in the system as a function of time (log scale) for all seven molecular systems at 200 K. b) Number of clusters per molecules in the system as a function of time (log scale) for all seven molecular at 200 K.In all simulations, the mechanism of aggregation is through small clusters coalescing into larger clusters. A possible alternative mechanism of a small seed which inducing single molecule aggregation was never observed. Two distinct growth phases can be observed as demonstrated in figure 2. In the

Figure 3 :Figure 4 :

 34 Figure 3: (a) Representation of naphthalene molecule with the vectors that define the φ angle and benzene molecule with the vectors that define the θ angle. Structural analysis of the orientation of molecules within a cluster for homogeneous systems for 200 K simulations as defined by two terms φ and θ as a function of distance, r. b-c) benzene-benzene interactions, d-e) naphthalene-naphthalene interactions and f-g) anthracene-anthracene interactions. The colours provide a measure of the number of molecules observed in the various configurations, from black indicating never observed to red indicating very frequent.

Figure 5 :

 5 Figure 5: Structural analysis of the orientation of molecules within a cluster for heterogenous systems at 200 K as defined by φ angle as a function of distance, r. Benzene-benzene, naphthalene-naphthalene and anthracene-anthracene interactions described in the presence of benzene, naphthalene and anthracene. The colours provide a measure of the number of molecules observed in the various configurations, with black indicating none and red indicating many interactions at that configuration with graduations from yellow to orange.

Figure 6 :

 6 Figure 6: Structural analysis of the orientation of molecules within a cluster for heterogenous systems as defined by φ angle as a function of distance, r, for benzene-benzene, benzene-naphthalene, benzene-anthracene, naphthalene-naphthalene naphthalene-anthracene and anthracene-anthracene interactions. The colours provide a measure of the number of molecules observed in the various configurations, with black indicating none and red indicating many interactions at that configuration with graduations from yellow to orange

Figure 7 :

 7 Figure 7: Final cluster of the three-species mixture for each of the five replica simulations. The centre of the cluster is predominantly made of anthracene molecules (orange), mostly surrounded by naphthalene (green) and benzene (grey) molecules in the outer surface.

Figure 8 .

 8 Figure 8. Parallel displaced (PD) (a) and T-shaped (b) configurations of dimers and trimers of benzene complexed with benzene, naphthalene, and anthracene. Equilibrium bond distances and PWPB95-D3/Def2-QZVPP binding energies are given (positive energies indicate an attractive interaction).

Table 1 :

 1 Composition of the mixture for the seven-system studied using molecular dynamics.

	System	n(Benzene) n(Naphthalene) n(Anthracene) n(atoms)
	Benzene	60	0	0	720
	Naphthalene	0	40	0	720
	Anthracene	0	0	30	720
	Benzene and Naphthalene	30	20	0	720
	Benzene and Anthracene	30	0	15	720
	Naphthalene and Anthracene	0	20	15	720
	Benzene, Naphthalene and	20	13	10	714
	Anthracene

  48,49

Acknowledgements

NAMD was developed by the Theoretical and Computational Biophysics Group in the Beckman Institute for Advanced Science and Technology at the University of Illinois at Urbana-Champaign. This work was in part supported by resources provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia. This research was undertaken in part with the assistance of computational resources from the Pople high-performance computing cluster of the Faculty of Science at the University of Western Australia and with the assistance of resources from the National Computational Infrastructure (NCI), which is supported by the Australian Government. AK gratefully acknowledges an Australian Research Council (ARC) Future Fellowship (Project No. FT170100373). I. S.-M. fellowship is funded by Australian Research Council (No. FT140100191). I.C. acknowledges funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC2013-SyG, Grant Agreement No. 610256 NANOCOSMOS. DT acknowledges MEAE and MESR for the travel support.

Conflict of Interest

The authors declare no conflict of interest.