Welding dynamics in an atomistic model of an amorphous polymer blend with polymer–polymer interface

We consider an atomistic model of thermal welding at the polymer-polymer interface of a polyetherimide/polycarbonate blend, motivated by applications to 3D manufacturing in space. We follow diffusion of semiflexible chains at the interface and analyze strengthening of the samples as a function of the welding time tw by simulating the strain-stress and shear viscosity curves. The time scales for initial wetting, and for fast and slow diffusion, are revealed. It is shown that each component of the polymer blend has its own characteristic time of slow diffusion at the interface. Analysis of strainstress demonstrates saturation of the Young’s modulus at tw = 240 ns, while the tensile strength continues to increase. The shear viscosity is found to have a very weak dependence on the welding time for tw > 60 ns. It is shown that both strain-stress and shear viscosity curves agree with experimental data.


INTRODUCTION
Understanding the properties of polymer-polymer interfaces represents a long-standing problem that is of both fundamental and technological importance [1][2][3] . For example, reptation, entanglement, and stress relaxation 1,4-8 determine the welding dynamics and parts strength in fused deposition modeling [9][10][11] , which is of significant interest for space applications including NASA's In-Space Manufacturing (ISM) project that seeks to develop the materials and processes needed to provide an on-demand manufacturing capability for deep space exploration missions 12 .
It is also known that the atomistic structure of polymers can substantially influence fundamental properties of welding dynamics mentioned above [13][14][15] . In particular, polyimides [15][16][17] and their blends [18][19][20][21] have recently attracted considerable attention due to the combination of processing, thermal, and mechanical properties suitable for aerospace applications. And the ISM project seeks to expand its on-orbit printing and recycling ecosystem by including blends of polyetherimide into the process. However, many issues related to optimization and predicting the quality of the parts remain unsolved [19][20][21][22] and their resolution could be facilitated by atomistic insight into the dynamics of polymer chains at an interface.
The latter simulations yield insight but do not, however, take into account limited flexibility and electrostatic properties of polymer chains that are of fundamental importance for chain diffusion 16,[31][32][33] . An additional challenge is posed by the need to understand the structure-properties relationship at the interface of polymer blends. Although MD has been widely used to investigate bulk properties of blends 25,[34][35][36] its application to modeling interfaces in blended polymer materials remains very limited. In particular, despite the fact that polyetherimide (PEI) polycarbonate (PC) mixtures are of significant importance in aerospace applications 37,38 , there have been only a few MD studies of their bulk and interfacial properties. For example, Zhang and co-authors 27  To address these issues we develop a fully atomistic model of the polymer-polymer interface in PEI/PC mixtures. The simulation of diffusing chains at the interface and relaxation of their distribution to the bulk state was limited to 240 ns. The welding analysis in this model has shown that initial wetting is followed by fast diffusion with Eyring's type jumps and then by slow diffusion of center-of-mass of polymer chains. The latter process was dominated by reptation with two characteristic timescales associated with two components of the blend. Existence of two slow timescales results in faster equilibration of the PC chains distribution, as compared to that of polyetherimide and in turn affects the dynamics of strengthening at welded interface.
The sample strength is characterized by simulation of uni-axial and shear deformations after quenching the samples to room temperature at t w = 60 and 240 ns. It is shown that the Young's modulus of the quenched samples increases between t w = 60 and t w = 240 ns and saturates at t w = 240 ns while, at the same time, the yield strength is continuing to increase. We further demonstrate that the dependence of the shear viscosity η on the shear rate in our model is linear on a log-log scale, corresponding to the expected shear-thinning behavior of the PEI/PC blends.
Both strain-stress curves and dependence of η on shear rate are shown to be in reasonable agreement with available experimental data. Thermal cycling was performed for an additional 1 µs to analyze the thermomechanical properties of welded samples, which were shown to be in agreement with experimental data, as will be discussed in detail elsewhere 39 .
The model size was limited to 61912 atoms, which in turn imposes limitations on the cell size and on the timescale of the analysis. In this sense, the results presented should be considered complementary to those obtained for large coarse-grained models.
The paper is organized as follows. In the next section we describe the model of polymerpolymer interface. The "Results" section discusses the interfacial diffusion of polymer chains, uni-axial sample deformation and the dependence of the shear viscosity on temperature and shear rate. Finally, the results are summarized and possible future work is outlined in the "Conclusions".

MODEL
According to Cicala et al. 38 the choices of high performance filaments for fused deposition modeling include Ultem 1000, Ultem 9085 and Polyphenylsulfone. Both Ultem grades are based on polyetherimide and they are certified for use in the automotive, medical and aerospace fields. Ultem 1000 is a pure PEI, while Ultem 9085 is a mixture of ∼80% PEI and ∼20% PC co-polymer blend incorporated for improved flow 20 . The present work focuses on the welding of materials such as Ultem 9085, studied through molecular dynamics simulations.
To model interface welding we prepared samples that consist of two amorphous cells with a mixture of polyetherimide and polycarbonate chains. Each cell had one flat face, see  The degree of polymerization (DP) at present is limited by the maximum size of the samples that could be processed in simulations. The resulting DP corresponds to the transition from oligomers to polymers. We note that the characteristic ratios for PC and PEI are 3-4 as estimated by Bicerano's group contribution method. The number of backbone bonds are 20 (PC) and 25 (PEI) for the models with the DP 5 meaning that the model polymers used in this study can behave like random polymer chains.
In the absence of experimental data demonstrating the dependence of the results on the DP we estimated this dependence using Fox-Flory theory 40 . According this theory both cases 5-5 and 6-8 (PEI-PC DP) are in the range ≥1000 g/mol where the dependence of the glass transition temperature T g on the DP is approaching a plateau. In addition, the MD computation 27 of the solubility parameter χ of PEI and PC chains as a function of the DP indicates that the change from 5-5 to 6-8 (PEI-PC DP) in the model corresponds to the transition towards saturation of the χ. Comparison of the two models reveals similar quantitative features of the interface diffusion and allows us to assume that the captured features are robust.
Each cell was built using the software package J-OCTA 41 by placing at random polycarbonate and polyetherimide chains such that initial dilute mixture had 20 wt % PC and 80 wt % PEI and the initial sample density ≈ 0.65 g/cm 3 . To prepare one nearly atomically flat surface (with normal vector parallel to the z-axis) in each cell we applied the Lennard-Jones (LJ) potential at this face while keeping the boundary conditions periodic in the X-and Y -directions and free in the Z-direction. The cells were relaxed compressed. and relaxed again using standard procedures 41 described in further details in section "Simulation details" The relaxation process was monitored by calculating potential energy, energy of nonbonding interaction, density, and radius of gyration of the PEI and PC chains as a function of time. These quantities level out by the end of the relaxation process. Note that in the actual manufacturing process the polymers at the interface are away from equilibrium (see more detailed discussion in the SI) and further analysis of the dependence of the interface diffusion on the distribution of the chains at the interface will be required in the future.
Finally, the interface welding was simulated by allowing polymer chains to diffuse freely at high temperature as shown in Fig. 1, see Table S4 of the SI. Top (cyan) and bottom (gray) cells in these sample were prepared as discussed above. We note that there is a second interface due to periodic boundary conditions in the merged system. However, initial entanglement at the second boundary is the same as in the bulk and initial stress relaxes fast to the bulk state. Welding takes place at the interface flat separated faces with initial gap shown in Fig. 1. and in what follows we focus on the analysis of this interface.

Simulation details
The welding simulations were conducted in an NPT ensemble with periodic boundary conditions using either LAMMPS 42,43 or GROMACS [44][45][46]  The limited size of the samples and periodic boundary conditions allow analysis of welding during a few hundred nanoseconds. In particular, simulations of welding in the 1st sample were performed at 600 K during 240 ns. After 60 ns and 240 ns, partially equilibrated samples were quenched to 300 K in 12 steps of 25K each. Additional thermal cycling of the smaller sample was performed between 300 K and 600 K with the time step varying between 12 ns and 25 ns. The initial temperature of the interface was chosen to correspond roughly to the actual welding process, in which an extrusion temperature is 623-653 K (i.e. well above T g = 459 K) while the bed temperature is ∼ 413 K (i.e. below T g ). 6 The temperature profile of the simulations is shown in Fig. S5 of the SI. The total time above the T g during thermal cycling was ∼400 ns. Welded samples at different temperatures obtained during thermal cycling were used to estimate thermomechanical and spectral properties of the PEI/ PC blends as discussed in a separate paper 39 .
The welding, quenching, and thermal cycling used an NPT ensemble i.e. keeping pressure, temperature, and the total number of particles fixed 1 . The quenched samples at time instances 60, 240 ns, and after thermal cycling with different thickness of the welded layer were used to study strain-stress and shear-viscosity -shear-rate relations of the blends.
The classical Hamiltonian of the polymer system used in this study is of the form 49 where the first term K on the right-hand side is the kinetic energy of the system, the next term U bond corresponds to bonding interactions between atoms in the polymer chain, and the last term corresponds to non-bonding interactions. Here r ij is the distance between the i-th and j-th atoms, q i is the atomic charge, ε 0 is the permittivity of free space, ε is the dielectric constant, and A ij and B ij are parameters of the LJ potential. In the simulations we use the Dreiding force field 50 and assign partial charges using the molecular orbital (MO) method with PM3 in MOPAC 51 , see the SI for further details. Electrostatic interactions are calculated by the particle-mesh technique, particle-mesh Ewald (PME) 52 for GROMACS and particle-particle-particle-mesh Ewald 53 for VSOP and LAMMPS.
We note that Dreiding force field was used earlier 27,54 for molecular dynamic simulations of the polyetherimide. In 27 its performance was compared to that of the COMPASS 55 force field. It was found that Dreiding gave the most accurate estimation of solubility parameters of the present systems and provided much faster performance for simulation work on polymer blends. It was also successfully used in 25,56 to simulated epoxy thermosets with related structure of the polymer chains. However, this force field underestimates density by 10-15 %, see Fig. S6 of the SI, and to use Dreiding force field in this work we performed extensive validation of other thermo-mechanical properties of PEI/PC blends and found a reasonable agreement with available experimental data 57 .
There are promising alternatives to Dreiding force field in terms of efficiency and accuracy 28,58,59 . The details of the validation and comparison with the performance of alternative force field will be reported elsewhere 39 .
To estimate the strengthening of samples as a function of welding time we note that in the presence of deformation the internal energy of the system E = ⟨H⟩ is related to the Gibbs G and Helmholtz F free energies as 60 where S is the entropy, T is the temperature, and σ ik and are respectively the stress and strain tensors given by 60 To account for the interface contribution to the overall stress we note also that the system consists of two bulk polymer samples and the interface. In general, the energy of such systems is a sum of three contributions 61 where the L, R, and s indices correspond respectively to the left, right, and surface components. In MD, the 6 components of the pressure tensor (negative of stress tensor) are calculated as 62 and f ik is the total force acting on the i-th atom.
Thus we see that the interface energy contributes to the value of the pressure tensor measured for the whole system. Changes of interfacial energy as a function of welding time can be detected by applying uni-axial or uni-diagonal deformation and measuring the corresponding pressure tensor. The gradient of the free energy at the sample interface is also a driving force for the diffusion that underlies the welding process. We will now discuss the welding dynamics observed in the simulations.

Diffusion at the interface
To resolve interfacial diffusion at the initial time we used the second, larger, sample with longer chains and with total number of atoms 61912. The interface was prepared and relaxed at a temperature 650 K, following the procedure described in the Model section. The interfacial dynamics in both samples was qualitatively similar. shows the initial dynamics of the non-binding energy in the large sample corresponding to transition from "wetting' to "slow" diffusion regime.
As the samples are allowed to equilibrate, the first phenomenon observed is "wetting", when the two surfaces quickly come close to each other 1,63 , on the time scale of a few picoseconds. The wetting process is governed by the electrostatic and van der Waals forces. The smaller the gap between two surfaces (assuming that Pauli repulsion for atoms in different chains remains weak), the smaller is the total non-bonding energy of the sample, as shown in the inset of Fig. 2, see also video of "wetting" process in the SI.
Next, we observe a relatively fast diffusion of polymer chains on a time scale of ∼20-30 ps. The difference in the diffusion rates between "slow" and "fast" diffusion can be clearly seen in the figure as the difference between two slopes shown by dashed lines. We attribute the observed accelerated diffusion to the initial existence of un-equilibrated chain ends and "vacancies" on the both sides of the interface. As a result the diffusion is driven by both reptation and Eyring-type 64 jumps of the chain ends between quasi-equilibrium positions: see the SI for further illustration of this point.  This extension of the distributions tails has a profound effect on the interface strength.
Indeed, according to Wool 6 full strength is obtained when polymer filaments are inter-diffused to a distance equal to 81% of the radius of gyration (R g ). (Note, that for the semiflexible chains studied here Wool's criterion should be used only as a guiding approximation, see further discussion in the SI.) This condition means that the centers-of-mass distribution of the chains at the interface should approach its bulk value. In the large sample, the PEI chains extend from their center of mass by ∼60Å, while the PC chains extend approximately ∼30Å, see Fig. S10 and S11 of the SI.
The overlapping of cores of the distributions by ∼81 % of the R g assumes that the maximum of the distribution for PEI chains approaches the interface by the distance ∼ 0.6R g .
Therefore, we expect complete healing of the interface when the maximum (relative to the interface located at z = 0Å) chain extension in the z-direction is ∼40 Å. We observe, however, a substantial slowing down of the tail extension beyond 20Å, which is attributable to the blend structure, with 80% of PEI and 20% of PC chains. The latter chains, being smaller and more flexible, diffuse faster towards the interface, while the stiffer and longer PEI chains need more time to equilibrate.
It is known that miscible 36   We conclude that the strength of the interface is approaching its bulk value as a function of time, but that the curing process remains incomplete. On the time scale of the simulations, the interface strength will be determined mainly by the inter-diffusion of PC chains and will be lower than that expected for polyetherimide. We now discuss the estimation of sample strength as a function of welding time, using MD simulations.

Strain-stress curve
The strain-stress curves of the samples were estimated by simulating uni-axial deformation at a constant rate in the Z-axis direction using a scenario developed by J-OCTA 41 . These simulations were performed in LAMMPS 65 using an NVT ensemble. The sample shape before and after the deformation is shown in Fig. S9 of the SI. The estimates are found to depend on the elongation rate v e . The overall conclusion of this section is that fully atomistic MD simulations yield esti- The following features may be noted in the figure.
First, we observe that the dependence of η on the shear rate exhibits characteristic shear- Both E a and η int increase as functions of welding time. However, the η int contribution to the overall viscosity of the sample is relatively small and it decreases with time due to the interface equilibration. As a result the observed shift of η is less than the measurement error.
The dependence of η on temperature is shown in the inset of Fig. 8 for shear rate γ ≈ 8 × 10 8 1/s. It is evident that the temperature dependence of η is well resolved and that it exhibits the expected trend of decreasing as the temperature rises. To facilitate comparison of the estimated η values with experimental data, we have extended MD simulations to deformation rates of ∼ 6.9 × 10 6 1/s using GROMACS on Amazon Web Services. The results of the extended simulations are compared with experimental data 66 in Fig. 9. We note that the MD predictions are in even better agreement with the experimental data obtained for Ultem 9085 at room temperature 71 .
These results demonstrate that atomistic simulations are capable of quantitative estimation of the shear viscosity in PEI/PC blends but are weakly sensitive to changes in the interface values of η.

CONCLUSIONS
We have developed an atomistic model of polyetherimide/polycarbonate blends with planar polymer-polymer interfaces. It takes explicit account of electrostatic interactions and the semi-flexible nature of the chains. We used molecular dynamics simulations of this model to analyze diffusion and strengthening as a function of welding time during the first 300 ns.
It was shown that welding occurs in a number of steps. The initial gap at the interface between the two pieces of polymer was closed on a time scale of a few pico-seconds in the well-known "wetting" process 1,63 .
During the second step we observed fast interfacial diffusion, attributable to the initial existence of un-equilibrated chain ends and "vacancies" on the both sides of the interface. Strengthening of the channels as a function of welding time was analyzed by simulating uni-axial elongation of quenched (to 300 K) samples after 60 and 240 ns of equilibration and after additional thermal cycling of the samples between 300 and 600 K. It was shown that Young's modulus E is increased when the welding time is changed from 60 ns to 240 ns.
Changes of E are much less pronounced for further increase of welding time during thermal cycling while the yield strength of the samples continues to increase.
The observed features were attributed to the slow CM diffusion of PEI and PC chains.
The increase of yield strength during thermal cycling corresponds to the slow equilibration of the CM distribution of PEI chains in accordance with Wool's criterion 6 .
We note that the breakup in all samples, except the breakages obtained after thermal cycling, was due to the pulling out of chains at the interface, which remained the weakest part of the system. After thermal cycling, breakup could occur at different locations indicating that almost the full strength of the interface had been recovered.
The shear viscosity η was estimated using MD simulations of uni-diagonal deformation of the samples quenched to 300 K after 60 and 240 ns of welding. It was shown that both samples exhibited the shear-thinning behavior characteristic of polymer melts 69 . The dependence of η on the shear rate was shown to be in good agreement with available experimental data.
However, we observed only weak (within the error of estimation) dependence of η on welding time in the simulations.
The thermal, mechanical, and optical properties obtained for welded samples also exhibit good agreement with available experimental data as will be discussed in detail in a separate paper 39 .
In conclusion, the results obtained demonstrate that fully atomistic models can be used to make realistic estimates of the parameters of welded polymer interfaces, and to anchor continuous models of polymer-based manufacturing processes. In particular, using obtained results we conjecture that the strengths of interfaces in polyetherimide blends (and parts produced using additive manufacturing) may be improved by reducing the molecular weight of PEI chains and by broadening their molar mass distribution. However, further more detailed research will be needed to verify this conjecture. The revealed features of the polymer dynamics at the interface are characteristic for semiflexible chains with partial charges and the results obtained could be useful for further development of the theory for such polymers.
It is also expected that the modeling approach developed in this work will help to elucidate specific features of materials and enhance physics-based characterization of polymer parts manufactured in space under micro-gravity conditions.
The main limitation of our results relates to the relatively small sizes of the samples and polymer chains. Fully atomistic simulation of the interface welding in larger models will be performed in future work. 20