Issue 
A&A
Volume 554, June 2013



Article Number  A4  
Number of page(s)  12  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201321325  
Published online  28 May 2013 
Static compression of porous dust aggregates
^{1} Department of Astronomical Science, School of Physical Sciences, Graduate University for Advanced Studies (SOKENDAI) Mitaka, 1818588 Tokyo, Japan
email: akimasa.kataoka@nao.ac.jp
^{2} National Astronomical Observatory of Japan, Mitaka, 1818588 Tokyo, Japan
^{3} Institute of Low Temperature Science, Hokkaido University, Kita, 0600819 Sapporo, Japan
^{4} Department of Physics, Nagoya University, Nagoya, 4648602 Aichi, Japan
^{5} Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Meguro, 1528551 Tokyo, Japan
^{6} Planetary Exploration Research Center, Chiba Institute of Technology, Narashino, 2750016 Chiba, Japan
Received: 20 February 2013
Accepted: 11 March 2013
Context. In protoplanetary disks, dust grains coagulate with each other and grow to form aggregates. While these aggregates are growing by coagulation, their filling factor φ decreases to φ ≪ 1; however, comets, the remnants of these early planetesimals, have φ ~ 0.1. Thus, static compression of porous dust aggregates is important in planetesimal formation. However, the static compressive strength has only been investigated for relatively highdensity aggregates (φ > 0.1).
Aims. We investigate and find the compressive strength of highly porous aggregates (φ ≪ 1).
Methods. We performed threedimensional Nbody simulations of aggregate compression with a particleparticle interaction model. We introduced a new method of static compression: the periodic boundary condition was adopted, and the boundaries move with low speed to get closer. The dust aggregate is compressed uniformly and isotropically by themselves over the periodic boundaries.
Results. We empirically derive a formula of the compressive strength of highly porous aggregates (φ ≪ 1). We check the validity of the compressive strength formula for wide ranges of numerical parameters, such as the size of initial aggregates, the boundary speed, the normal damping force, and material. We also compare our results to the previous studies of static compression in the relatively highdensity region (φ > 0.1) and confirm that our results consistently connect to those in the highdensity region. The compressive strength formula is also derived analytically.
Key words: planets and satellites: formation / methods: numerical / methods: analytical / protoplanetary disks
© ESO, 2013
1. Introduction
Planetesimal formation is a key issue in the study of how planets form in protoplanetary disks (Hayashi et al. 1985; Weidenschilling & Cuzzi 1993). However, the collisional growth of the dust from submicronsized dust to kilometersized planetesimals is still unknown.
In the growth process, one of the most important but unresolved problems is to determine the internal structural evolution of dust aggregates. The internal structure of dust is important in planetesimal formation because the dynamics of dust aggregates in protoplanetary disks are determined by coupling of gas and dust, in other words, the size and internal density of dust aggregates. In the early stage of dust coagulation in protoplanetary disks, the collision energy of the aggregates is too low to cause collisional compression (Blum 2004; Ormel et al. 2007; Zsom et al. 2010, 2011; Okuzumi et al. 2012). As a result, the internal mass density ρ decreases to ρ < 1.0 g cm^{3}.
Both theoretical and experimental studies have shown that mutual collisions lead dust aggregates to have their fractal dimension D ~ 2, which is socalled ballistic clustercluster aggregation (BCCA; Smirnov 1990; Meakin 1991; Kempf et al. 1999; Blum & Wurm 2000; Krause & Blum 2004; Paszun & Dominik 2006). The dust aggregates would be gradually compacted or disrupted in coagulation because of the increase in impact energy. This compaction has been investigated with numerical Nbody simulations that consider particleparticle interactions (Dominik & Tielens 1997; Wada et al. 2007, 2008, 2009; Suyama et al. 2008, 2012; Paszun & Dominik 2008, 2009; Seizinger et al. 2012).
Material parameters in our simulation.
In most previous studies investigating dust growth in protoplanetary disks, dust grains have been assumed to have constant internal mass density for simplicity (Nakagawa et al. 1981; Tanaka et al. 2005; Brauer et al. 2008; Birnstiel et al. 2010). However, dust porosity evolves during dust growth in real protoplanetary disks. In recent dust coagulation calculations, porosity evolution has been considered to be based on experimental and theoretical results (Ormel et al. 2007; Okuzumi et al. 2009, 2012; Zsom et al. 2011). These results also suggest that ρ decreases as ρ ≪ 0.1 g cm^{3}.
In the most recent work, the dominant coagulation mode has been shown to be similarsize collisions of dust aggregates though dust grains have size distribution (Okuzumi et al. 2012). As a result, their fractal dimension is approximately equal to two, and their internal mass density ρ has been shown to become 10^{5} g cm^{3} (equivalent to be the filling factor φ = 10^{5} for ice particles with a density of 1.0 g cm^{3}). Such fluffy dust aggregates are believed to become planetesimals. Since comets in our solar system, which would be remnants of planetesimals, have their internal mass density of ~0.1 g cm^{3} (A’Hearn 2011), dust aggregates must be compressed from ρ ≪ 0.1 g cm^{3} to ρ ~ 0.1 g cm^{3} in protoplanetary disks.
Compression at dust aggregate collisions has been investigated in previous studies. When collisional impact energy exceeds the critical energy, dust aggregates are compacted by their collision (e.g., Dominik & Tielens 1997; Suyama et al. 2008; Wada et al. 2007, 2008, 2009). However, the collisional compression is not effective at compressing dust aggregates (Okuzumi et al. 2012).
One of the other compression mechanisms in protoplanetary disks is static compression by disk gas or selfgravity. The static compressive strength of dust aggregates has been investigated both experimentally and numerically (Paszun & Dominik 2008; Güttler et al. 2009; Seizinger et al. 2012). However, the compressive strength has examined only relatively compact aggregates with ρ ≳ 0.1 g cm^{3} because their initial aggregates are ballistic particlecluster aggregation (BPCA) clusters. Because ρ decreases to ρ ≪ 0.1 g cm^{3}, at least in the early stage of dust growth, we need to reveal the static compressive strength with ρ ≪ 0.1 g cm^{3}.
In this work, we investigate the static compression of highly porous aggregates with ρ < 0.1 g cm^{3} by means of numerical simulations and an analytical approach. It is challenging to perform numerical simulations of the static and uniform compression of highly porous aggregates. Because such porous aggregates have low sound speed, we have to compress them at a much slower velocity than in the case of compact aggregates, as is shown in our simulations. Such a slow compression of the fluffy aggregates costs much computational time.
In previous numerical studies of static compression, a dust aggregate is compressed by a wall moving in one direction (Paszun & Dominik 2008; Seizinger et al. 2012). However, this method has disadvantages when reproducing uniform and isotropic compression. There are also side walls that do not move. These side walls also obstruct the tangential motion of monomers in contact with the walls, causing artificial stress on the aggregate, which restructures them. Moreover, since they measure the pressure with the force on the moving wall, the side walls may affect the pressure measurement. In the present work, we develop a new method reproducing static compression. Instead of the walls, we adopt periodic boundary conditions and the boundaries get closer to each other. With these slowly moving periodic boundaries, the aggregate is compressed uniformly and naturally. The periodic boundary condition also enables us to represent a much larger aggregate than inside the computational region. This saves on computational time remarkably.
This paper is organized as follows. We describe the model of our numerical simulations in Sect. 2. We show the results of our simulations and find the compressive strength in Sect. 3. We confirm the obtained compressive strength formula analytically in Sect. 4, and present our conclusion in Sect. 5.
2. Simulation setting
We performed threedimensional numerical simulations of the compression of a dust aggregate consisting of a number of spherical monomers. As the initial aggregate, we adopted a BCCA cluster. In this method, we solved interactions between all monomers in contact in each time step. Interactions between monomers in contact are formulated by Dominik & Tielens (1997) and reformulated by using the potential energies by Wada et al. (2007). We used the interaction model proposed by Wada et al. (2007) in this work. We briefly summarize the particle interaction model and material constants (see Wada et al. 2007, for details). Moreover, we describe the additional damping force in normal direction and the simulation setting in this section. In our simulations, the aggregate is gradually compressed by its copies over the moving periodic boundaries. This is an appropriate method of simulating uniform and isotropic compression. We also describe the boundary condition in this section. Since we do not have walls to measure the pressure in the periodic boundary condition, we use a similar manner of pressure measurement in molecular dynamics simulations. We also introduce the method of pressure measurement below.
2.1. Interaction model
We calculate the direct interaction of each connection of particles, taking all mechanical interactions modeled by Dominik & Tielens (1997) and Wada et al. (2007) into account. The material parameters are the monomer radius r_{0}, surface energy γ, Young’s modulus E, Poisson’s ratio ν, and the material density ρ_{0}. Table 1 lists the values of the material parameters for ice and silicate.
We perform Nbody simulations with ice particles except for one case with silicate particles. In protoplanetary disks, ice particles are the most dominant dust material beyond the snowline. Moreover, the computational time required for calculating ice particles is less than for silicate. Thus, we adopt ice particles in most simulations. We also treat a silicate case to compare with a previous study (Seizinger et al. 2012).
The critical displacement still shows a discrepancy between theoretical (ξ_{crit} = 2 Å) and experimental (ξ_{crit} = 32 Å) studies (Dominik & Tielens 1997; Heim et al. 1999). We adopt the same parameter as in Wada et al. (2011), ξ_{crit} = 8 Å as a typical length for ice particles, and ξ_{crit} = 20 Å for silicate particles to compare with Seizinger et al. (2012).
The parameter ξ_{crit} is related to strength of rolling motion. The rolling motion between monomers is crucial in compression. The rolling energy E_{roll} is the energy required to rotate a particle around a connecting point by 90°. The rolling energy can be written as $\begin{array}{ccc}{\mathit{E}}_{\mathrm{roll}}\mathrm{=}\mathrm{6}{\mathit{\pi}}^{\mathrm{2}}\mathit{\gamma}{\mathit{r}}_{\mathrm{0}}{\mathit{\xi}}_{\mathrm{crit}}& & \end{array}$(1)(see Wada et al. 2007, for details). In the case of ice monomers, for example, E_{roll} = 4.37 × 10^{9} erg for ξ_{crit} = 8 Å.
We use a normalized unit of time in our simulations. For ice particles, the normalized unit of time is $\begin{array}{ccc}{\mathit{t}}_{\mathrm{0}}\mathrm{=}\mathrm{0.95}\left(\frac{{\mathit{\rho}}_{\mathrm{0}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{r}}_{\mathrm{0}}^{\mathrm{7}\mathit{/}\mathrm{6}}}{{\mathit{E}}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{\gamma}}^{\mathrm{1}\mathit{/}\mathrm{6}}}\right)\mathrm{=}\mathrm{1.93}\mathrm{\times}{\mathrm{10}}^{10}\mathrm{s}\mathit{,}& & \end{array}$(2)which is a characteristic time, and it approximately represents the oscillation time of particles in contact at the critical collision velocity (see Wada et al. 2007, for details).
2.2. Damping force in normal direction
The normal force between two monomers is repulsive when the monomers are close or attractive when they are stretched out. Thus, normal oscillations occur at each connection. For realistic particles, these oscillations would dissipate because of viscoelasticity or hysteresis in the normal force (e.g., Greenwood & Johnson 2006; Tanaka et al. 2012). For such damping of normal oscillation, we add an artificial normal damping force to the particle interaction model, following the previous studies (Suyama et al. 2008; Paszun & Dominik 2008; Seizinger et al. 2012).
Assuming that two particles in contact have their position vectors x_{1} and x_{2}, respectively, the contact unit vector n_{c} is defined as $\begin{array}{ccc}{{n}}_{\mathrm{c}}\mathrm{=}\frac{{{x}}_{\mathrm{1}}\mathrm{}{{x}}_{\mathrm{2}}}{\mathrm{}{{x}}_{\mathrm{1}}\mathrm{}{{x}}_{\mathrm{2}}\mathrm{}}& & \end{array}$(3)(see Fig. 2 in Wada et al. 2007). We introduce a damping force between contact particles in normal direction, defined as $\begin{array}{ccc}{{F}}_{\mathrm{damp}}\mathrm{=}\mathrm{}{\mathit{k}}_{\mathrm{n}}\frac{{\mathit{m}}_{\mathrm{0}}}{{\mathit{t}}_{\mathrm{0}}}{{n}}_{\mathrm{c}}\mathrm{\xb7}{{v}}_{\mathrm{r}}\mathit{,}& & \end{array}$(4)where k_{n} is the damping coefficient in normal direction and m_{0} the monomer mass. The adopted value of k_{n} is on the order of 0.01. To show that the result is independent of the normal oscillation damping, we perform Nbody simulations with the damping factor k_{n} as a parameter.
The timescale of damping is $\begin{array}{ccc}{\mathit{\tau}}_{\mathrm{damp}}\mathrm{~}\frac{{\mathit{t}}_{\mathrm{0}}}{{\mathit{k}}_{\mathrm{n}}}\mathrm{~}{\mathrm{10}}^{\mathrm{2}}{\mathit{t}}_{\mathrm{0}}\mathit{,}& & \end{array}$(5)for k_{n} = 0.01, it is much shorter than the simulation timescale, which is typically ~10^{7}t_{0}. We show that the obtained compressive strength is independent of the artificial normal damping force in our simulations (see Sect. 3.4).
2.3. Uniform compression by moving boundaries
We adopt the periodic boundary condition in our simulations. The aggregate in the computational region is surrounded by its copies, as shown in Fig. 1.
Fig. 1 Schematic drawing of the periodic boundary condition. Each box illustrates a boundary box with a side length L for all directions. When the boundary starts to get closer, the aggregate sticks to the neighboring aggregates over the boundary and is compressed by them. It should be noted that this picture is illustrated in the 2D direction, but our simulations are performed in 3D. 
Initially, we set a cubic box whose sides are periodic boundaries with a size of L larger than the aggregate. Thus, the initial BCCA cluster is detached from its neighboring copies over the periodic boundaries. In our simulations, we gradually move the boundaries to the center of the aggregate to get closer to one another. As a result, the aggregate sticks to the neighboring copies and is compressed by them in a natural way. Therefore, the aggregate in the computational region corresponds to a small part of a whole large aggregate. In other words, although the number of particles in numerical simulations are limited because of computational cost, the periodic boundary condition enables us to investigate a large aggregate, such as a ~ cmsized dust aggregate in protoplanetary disks.
Another advantage of the periodic boundary condition is that we do not need to introduce the wall for compression. In the previous Nbody simulations of static compression, dust aggregates are compressed by using the wall against the dust aggregate (Paszun & Dominik 2008; Seizinger et al. 2012). The wall itself may have some artificial effects on such experiments. For example, the wall moves in one direction and thus this may be different from isotropic compression. Besides, wallparticle interaction is different from particleparticle interaction, so it must be treated carefully. In contrast, the periodic boundary condition does not need walls for compression because a dust aggregate is compressed by the neighboring aggregate over the periodic boundary. In addition, the periodic boundaries in three directions make it possible to compress the aggregate isotropically. We calculate not only the interactions of particles in contact inside the computational region but also the interactions of the particles in contact across the periodic boundaries. Thus, no special treatment of interactions, which is wallparticle interactions in the case of simulations with walls in previous studies, is required when a particle crosses the periodic boundaries.
The computational cubic region has length L, and the coordinates in x,y, and z directions are set to be − L/2 < x < L/2, − L/2 < y < L/2, and − L/2 < z < L/2, respectively. We adopt periodic boundary conditions for all directions to reproduce a part of a large aggregate, where L decreases with time t, L = L(t). The initial size of the box L_{0} is adopted as the maximum size of the dust aggregate in x, y, and z directions.
With the settings above, we move the boundaries of the computational region toward the center of the region. The velocity at the boundary is given by $\begin{array}{ccc}{\mathit{v}}_{\mathrm{b}}\mathrm{=}\mathrm{}\frac{{\mathit{C}}_{\mathrm{v}}}{{\mathit{t}}_{\mathrm{0}}}\mathit{L}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}& & \end{array}$(6)where C_{v} is a dimensionless parameter (we call C_{v} the strain rate parameter hereafter). Owing to this definition of the boundary speed, the aggregate is compressed at a constant strain rate independent of the region scale L.
The box size decreases with the constant rate C_{v} in three directions. This corresponds to isotropic compression. Since $\frac{\mathrm{d}\mathit{L}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{2}{\mathit{v}}_{\mathrm{b}}$, the box size is written as $\begin{array}{ccc}\mathit{L}\mathrm{=}{\mathit{L}}_{\mathrm{0}}\mathrm{exp}\left(\mathrm{}\mathrm{2}{\mathit{C}}_{\mathrm{v}}\frac{\mathit{t}}{{\mathit{t}}_{\mathrm{0}}}\right)\mathrm{\xb7}& & \end{array}$(7)Therefore, the whole time of compression is t_{0}/C_{v}. Typically we chose C_{v} = 3 × 10^{7} so the compression time is ~3 × 10^{6} t_{0} ~ 0.6 ms.
When a particle crosses a periodic boundary, the velocity should be treated carefully to reproduce the quasistatic compression with periodic boundary condition. Figure 2 illustrates how to calculate the velocity of particles across the periodic boundary.
Fig. 2 Schematic drawing to illustrate how the particle velocity is calculated when a particle crosses a periodic boundary. For simplicity, we consider this situation in a 2D field, but we actually calculate this in a 3D situation. We consider that a dust particle is close to the boundary in the left figure. In the next time step, the particle crosses the boundary (dashed circle in the right figure). We put the particle on the other side of the boundary as expressed in Eqs. (8) and (10). The velocity component is converted as expressed in Eqs. (9) and (11). This treatment reproduces the isotropic compression in the velocity field well. 
When a particle goes out of the computational region across the boundary at x = L/2, we relocate the particle to the opposite side (i.e., from the boundary at x = −L/2). In that case, the position of the particle in x direction is converted as $\begin{array}{ccc}\mathit{x}\mathrm{}\mathrm{}\mathrm{\to}\mathit{x}\mathrm{}\mathit{L}& & \end{array}$(8)Since the two boundaries at x = − L/2 and x = L/2 have a relative velocity of 2v_{b}, the xcomponent of the velocity v_{x} of the particle is also converted as $\begin{array}{ccc}{\mathit{v}}_{\mathit{x}}\mathrm{}\mathrm{}\mathrm{\to}{\mathit{v}}_{\mathit{x}}\mathrm{+}\mathrm{2}{\mathit{v}}_{\mathrm{b}}\mathit{.}& & \end{array}$(9)Owing to the conversion of v_{x}, the velocity of particle against the boundary that the particle crosses does not change before and after the crossing. For a particle across the boundary at x = − L/2, the position and the velocity are converted as $\begin{array}{ccc}& & \\ & & \end{array}$\arraycolsep1.75ptWe also have the same treatments for particles across the boundaries at y = ± L/2 and z = ± L/2.
We introduce the constant strain rate at the boundaries for scaleless discussion. However, the initial aggregate is not moving. As the simulation starts, if all the particles in the aggregate are not moving, only the particles close to the boundaries have initial velocity. This is not a constant strain rate. To reproduce the scaleless constant strain rate initially, therefore, we first give all monomers the velocity smoothly connected to the boundary speed. The initial velocity is expressed as $\begin{array}{ccc}{v}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathrm{b}}\mathrm{\times}\frac{{r}}{{\mathit{L}}_{\mathrm{0}}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(12)where r is the position vector of the monomers.
Fig. 3 Snapshots of the evolution of an aggregate under compression in the case of N = 16 384. The top three figures are 3D visualizations. They have the same scale with different time epochs. The white particles are inside a box enclosed by the periodic boundaries. The yellow particles are in neighboring boxes to the box of white particles. For visualization, we do not draw the copies on the back and front sides of the boundaries but only 8 copies of the white particles across the boundaries. Each bottom figure represents projected positions onto 2D plane of all particles in each corresponding top figure. The gray points in the bottom figures correspond to the positions of the white particles in the top figures, and the yellow points correspond to those of the yellow particles in the top figures. Scales are in μm. 
2.4. Pressure measurement
In previous studies, a dust aggregate is enclosed by walls, and the pressure is calculated by measuring the force exerted on the walls by the dust aggregate. In this work, a dust aggregate is compressed by themselves because of the periodic boundary condition. Therefore, we introduce another method of measuring the pressure on the aggregate. We calculate the pressure of the dust aggregate in the standard way in molecular dynamics simulations using the virial theorem as follows (e.g., Haile 1992).
Fig. 4 a) Pressure P in [Pa] against filling factor φ. The ten thin solid lines show the results for the initial BCCA clusters with different initial random numbers and thick solid line shows the arithmetic average of the ten runs. b) Pressure P in [Pa] against filling factor φ. Same as the thick solid line in a) plotted with a dotted line of Eq. (25). The parameters are N = 16 384, C_{v} = 3 × 10^{7}, k_{n} = 0.01, and ξ_{crit} = 8 Å. 
We consider a virtual box that encloses the aggregate under consideration. We define the force acting from the walls of the virtual box on the particle i as W_{i}, and the sum of the forces from other particles on the particle i as F_{i}. The equation of motion of the particle i is given by $\begin{array}{ccc}\mathit{m}\frac{{\mathrm{d}}^{\mathrm{2}}{{r}}_{\mathit{i}}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\mathrm{=}{{W}}_{\mathit{i}}\mathrm{+}{{F}}_{\mathit{i}}\mathit{.}& & \end{array}$(13)We take a scalar product of both sides of the equation with r_{i} and take a long time average of both sides with time interval τ. The lefthand side becomes $\begin{array}{ccc}\mathit{m}\frac{\mathrm{1}}{\mathit{\tau}}{\mathrm{\int}}_{\mathrm{0}}^{\mathit{\tau}}{{r}}_{\mathit{i}}\mathrm{\xb7}\frac{{\mathrm{d}}^{\mathrm{2}}{{r}}_{\mathit{i}}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\mathrm{=}\mathit{m}\frac{\mathrm{1}}{\mathit{\tau}}\left[{{r}}_{\mathit{i}}\mathrm{\xb7}\frac{\mathrm{d}{{r}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}\right]\begin{array}{}\mathit{\tau}\\ \mathrm{0}\end{array}\mathrm{}\mathit{m}\frac{\mathrm{1}}{\mathit{\tau}}{\mathrm{\int}}_{\mathrm{0}}^{\mathit{\tau}}\frac{\mathrm{d}{{r}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}\mathrm{\xb7}\frac{\mathrm{d}{{r}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}\mathrm{d}\mathit{t}\mathit{.}& & \end{array}$(14)The first term on the righthand side vanishes in the limit of τ → ∞. We define the takingalongtime average in t as ⟨ ⟩ _{t}. Taking a summation of all particles and a long time average of Eq. (13), we obtain $\begin{array}{ccc}{\u27e8\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{\mathrm{1}}{\mathrm{2}}\mathit{m}{\left(\frac{\mathrm{d}{{r}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}\right)}^{\mathrm{2}}\u27e9}_{\mathit{t}}\mathrm{=}\mathrm{}\frac{\mathrm{1}}{\mathrm{2}}{\u27e8\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{{r}}_{\mathit{i}}\mathrm{\xb7}\mathrm{(}{{W}}_{\mathit{i}}\mathrm{+}{{F}}_{\mathit{i}}\mathrm{)}\u27e9}_{\mathit{t}}\mathit{.}& & \end{array}$(15)The first term on the righthand side is related to pressure P. The pressure is an average of all forces acting on the wall from all particles. Using the normal vector n of the wall surface directed outward, the force received by the wall that has an area dS is PndS. Therefore, $\begin{array}{ccc}{\u27e8\sum _{\mathit{i}}{{r}}_{\mathit{i}}\mathrm{\xb7}{{W}}_{\mathit{i}}\u27e9}_{\mathit{t}}\mathrm{=}\mathrm{}{\mathrm{\int}}_{\mathit{S}}\mathit{P}{n}\mathrm{\xb7}{r}\mathrm{d}\mathit{S}\mathrm{=}\mathrm{}\mathrm{3}\mathit{PV}\mathit{.}& & \end{array}$(16)This equation is obtained by taking surface integral as $\begin{array}{ccc}{\mathrm{\int}}_{\mathit{S}}{n}\mathrm{\xb7}{r}\mathrm{d}\mathit{S}\mathrm{=}{\mathrm{\int}}_{\mathit{V}}\mathrm{div}{r}\mathrm{d}\mathit{V}\mathrm{=}{\mathrm{\int}}_{\mathit{V}}\left(\frac{\mathit{\partial x}}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial y}}{\mathit{\partial y}}\mathrm{+}\frac{\mathit{\partial z}}{\mathit{\partial z}}\right)\mathrm{d}\mathit{V}\mathrm{=}\mathrm{3}\mathit{V}\mathit{.}& & \end{array}$(17)The translational kinetic energy K, averaged over a long time, is given by $\begin{array}{ccc}\mathit{K}\mathrm{=}{\u27e8\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{\mathrm{1}}{\mathrm{2}}\mathit{m}{\left(\frac{\mathrm{d}{{r}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}\right)}^{\mathrm{2}}\u27e9}_{\mathit{t}}\mathrm{\xb7}& & \end{array}$(18)Using K and P, Eq. (15) gives an expression of P as $\begin{array}{ccc}\mathit{P}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\mathit{K}\mathit{/}\mathit{V}\mathrm{+}\frac{\mathrm{1}}{\mathrm{3}}{\u27e8\sum _{\mathit{i}}{{r}}_{\mathit{i}}\mathrm{\xb7}{{F}}_{\mathit{i}}\u27e9}_{\mathit{t}}\mathit{/}\mathit{V}\mathit{.}& & \end{array}$(19)We define the force from particle j on particle i as f_{i,j}. Force F_{i} can be written as a summation of the force from another particle as $\begin{array}{ccc}{{F}}_{\mathit{i}}\mathrm{=}\sum _{\mathit{j}\ne \mathit{i}}{{f}}_{\mathit{i,j}}\mathit{.}& & \end{array}$(20)Using f_{i,j} = − f_{j,i}, we finally obtain the pressure measuring formula as $\begin{array}{ccc}\mathit{P}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\mathit{K}\mathit{/}\mathit{V}\mathrm{+}\frac{\mathrm{1}}{\mathrm{3}}{\u27e8\sum _{\mathit{i}\mathit{<}\mathit{j}}\mathrm{(}{{r}}_{\mathit{i}}\mathrm{}{{r}}_{\mathit{j}}\mathrm{)}\mathrm{\xb7}{{f}}_{\mathit{i,j}}\u27e9}_{\mathit{t}}\mathit{/}\mathit{V}\mathit{.}& & \end{array}$(21)The first term on the righthand side of the equation represents the translational kinetic energy per unit volume, and the second term represents the summation of the force acting at all connections per unit volume. This expression is useful for measuring the pressure of a dust aggregate under compression. We do not need to put any artificial object, such as walls, in simulations because Eq. (21) is totally expressed in terms of the summation of the physical quantities of each particle, which are the mass, the position, the velocity, and the force acting on the particle. In our calculations, we take an average of pressure for every 10 000 time steps, corresponding to 1000 t_{0} because we set 0.1 t_{0} as one time step in our simulation.
As mentioned in Sect. 2.2, the adopted damping force corresponds to rapid damping of normal oscillations. Thus, the kinetic energy of random motion rapidly dissipates. This corresponds to the static compression, and thus the compressive strength is determined by the second term of Eq. (21).
3. Results
The top three panels of Fig. 3 show snapshots of the evolution of an aggregate under compression in the case where N = 16 384, C_{v} = 3 × 10^{7}, k_{n} = 0, and ξ_{crit} = 8 Å. The top three panels have the same scale but different time epochs, which are t = 0, 1 × 10^{6}t_{0}, and 2 × 10^{6}t_{0}. The white particles are inside the computational region enclosed by the periodic boundaries, while the yellow particles are in the neighboring copy regions (for visualization, we do not draw particles on the front and backsides copy regions). The bottom three panels represent the projected positions onto the twodimensional plane for the correspondent top three figures. We confirm that the dust aggregate is compressed by their copies from all directions. As the compression proceeds, the aggregate of white particles is compressed by the neighboring aggregate of yellow particles. We focus on how high pressure is generated by quasistatic compression in numerical simulations. Our numerical simulations have several parameters: the size of the initial BCCA cluster, the compression rate, the normal damping force, and the critical displacement (corresponds to the rolling energy). We investigate the dependence of the pressure on these parameters, by performing several runs with different parameter sets. Although we assume ice aggregates in most runs, we also investigate cases of silicate aggregates to compare them with previous studies.
3.1. Fiducial run: obtaining the compressive strength
We put a BCCA cluster as the initial aggregate. For each set of parameters, we randomly create ten BCCA clusters following Okuzumi et al. (2009) and take arithmetic averages of the ten simulations of the different initial clusters. The pressure is measured using Eq. (21) at each run. We define the filling factor of an aggregate as $\begin{array}{ccc}\mathit{\phi}\mathrm{=}\frac{{\mathit{V}}_{\mathrm{0}}\mathit{N}}{\mathit{V}}\mathit{,}& & \end{array}$(22)where V_{0} is the monomer volume, N the number of monomers of the aggregate, and V the volume enclosed by the boundaries, which has a length of L. The filling factor also can be written as φ = ρ/ρ_{0}. Figure 4 shows that the measured pressure as a function of the filling factor φ(t). The parameters of the simulations are N = 16 384, C_{v} = 3 × 10^{7}, k_{n} = 0.01, and ξ_{crit} = 8 Å. The corresponding E_{roll} is 4.74 × 10^{9} erg for ξ_{crit} = 8 Å. Each colored line in Fig. 4a shows each simulation with the different initial shape of the aggregate. Figure 4b shows the arithmetic average of the pressure measured in ten different runs. Each line shows in different ranges of φ. The lowest φ is determined with the largest size of the initial boundary boxes of the ten runs. We find that the compressive strength is reproduced well by $\begin{array}{ccc}\mathit{P}\mathrm{=}{\mathit{P}}_{\mathrm{0}}{\mathit{\phi}}^{\mathrm{3}}\mathit{,}& & \end{array}$(23)where P_{0} = 4.74 × 10^{5} Pa. We analytically discuss why the compressive strength is proportional to φ^{3} in Sect. 4. In the highdensity region (φ ≳ 10^{1}), the measured strength deviates from the line of P = P_{0}φ^{3}. This is because the dissipation mechanism changes in the highdensity region (see Sect. 3.4). The deviation in the lowdensity region (φ ≲ 3 × 10^{3}) is partly caused by a finite boundary speed (or compression rate) as discussed in the next section. Another reason for the deviation in the lowdensity region is related to the density of the initial BCCA cluster. The filling factor of BCCA φ_{BCCA} is estimated as $\begin{array}{ccc}{\mathit{\phi}}_{\mathrm{BCCA}}\mathrm{=}\frac{{\mathit{V}}_{\mathrm{0}}\mathit{N}}{{\mathit{V}}_{\mathrm{BCCA}}}\mathrm{=}{\left(\frac{\mathrm{3}}{\mathrm{5}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{N}}^{\mathrm{}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(24)where we use the radius and the volume of a BCCA cluster, ${\mathit{r}}_{\mathrm{BCCA}}\mathrm{=}\sqrt{\mathrm{5}\mathit{/}\mathrm{3}}{\mathit{N}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{r}}_{\mathrm{0}}$ and ${\mathit{V}}_{\mathrm{BCCA}}\mathrm{=}\mathrm{(}\mathrm{4}\mathit{\pi}\mathit{/}\mathrm{3}\mathrm{)}{\mathit{r}}_{\mathrm{BCCA}}^{\mathrm{3}}$, respectively (e.g., Suyama et al. 2008). For N = 16 384, we obtain φ_{BCCA} ~ 3 × 10^{3}. In the early stage of compression, φ is lower than φ_{BCCA} because the initial BCCA clusters are apart from each other. This space between BCCA clusters would also cause the deviation from the line of P = P_{0}φ^{3}.
We now discuss the coefficient P_{0} of the compressive strength. Wada et al. (2008) show that E_{roll} is important in the collisional compressive strength. Thus, E_{roll} is expected to also be important in the static compressive strength. Considering that the characteristic volume is the monomer’s volume ~ ${\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}$, we suppose ${\mathit{P}}_{\mathrm{0}}\mathrm{=}{\mathit{E}}_{\mathrm{roll}}\mathit{/}{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}$, based on dimension analysis. Therefore, the compressive strength can be written as $\begin{array}{ccc}\mathit{P}\mathrm{=}\frac{{\mathit{E}}_{\mathrm{roll}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}}{\mathit{\phi}}^{\mathrm{3}}\mathrm{\xb7}& & \end{array}$(25)We analytically discuss and confirm this equation in Sect. 4. We also plot this equation in Fig. 4b. This figure clearly shows that the result is well fit by Eq. (25).
We show that compressive strength is proportional to ξ_{crit}, which is proportional to the rolling energy E_{roll} in Sect. 3.5. We also confirm that Eq. (25) is applicable to the case of different r_{0} in the silicate case.
3.2. Dependence on the boundary speed
To statically compress the aggregate, we should move the boundary at a low enough velocity not to create inhomogeneous structure. Figure 5 shows the dependency on the strain rate parameter.
Fig. 5 Pressure P in [Pa] against filling factor φ with different strain rate parameter C_{v}. Each line shows the average of ten runs of the fixed strain rate: C_{v} = 1 × 10^{7}, 3 × 10^{7}, 1 × 10^{6}, 3 × 10^{6}, 1 × 10^{5}. The other parameters are the same for every ten runs: N = 16 384, k_{n} = 0.01, and ξ_{crit} = 8 Å. The dashed line is Eq. (25). 
Each line shows the average of ten runs. The fixed parameters are N = 16 384, k_{n} = 0.01, and ξ_{crit} = 8 Å. The strain rate parameter C_{v} is equal to 1 × 10^{7}, 3 × 10^{7}, 1 × 10^{6}, 3 × 10^{6}, and 1 × 10^{5}. The higher C_{v}, the higher pressure in the lowdensity region is required for compression. This is mainly caused by the ram pressure from the boundaries with high speed.
When the compression proceeds and the density becomes higher to reach the line of Eq. (25), the pressure follows the equation. From Fig. 5, C_{v} = 3 × 10^{7} creates a sufficiently low boundary speed. The boundary speed can be calculated as a function of φ. Using Eq. (6) and $\mathit{\phi}\mathrm{=}\mathrm{(}\mathrm{4}\mathit{/}\mathrm{3}\mathrm{)}\mathit{\pi}{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}\mathit{N}\mathit{/}{\mathit{L}}^{\mathrm{3}}$, the velocity difference between a boundary and the next boundary, v_{d}, can be written as $\begin{array}{ccc}{\mathit{v}}_{\mathrm{d}}\mathrm{=}\mathrm{\left}\mathrm{2}{\mathit{v}}_{\mathrm{b}}\mathrm{\right}\mathrm{=}\mathrm{2}\frac{{\mathit{C}}_{\mathrm{v}}}{{\mathit{t}}_{\mathrm{0}}}{\left(\frac{\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi}{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}\mathit{N}}{\mathit{\phi}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{\xb7}& & \end{array}$(26)In the case of C_{v} = 3 × 10^{7}, v_{d} = 12.7, 5.9, and 2.7 cm/s for φ = 10^{3}, 10^{2}, and 10^{1}, respectively.
Here, we discuss the velocity difference of boundaries, comparing with the effective sound speed of the aggregates. The effective sound speed can be estimated as $\begin{array}{ccc}{\mathit{c}}_{\mathrm{s}\mathit{,}\mathrm{eff}}\mathrm{~}\sqrt{\frac{\mathit{P}}{\mathit{\rho}}}\mathrm{~}\sqrt{\frac{{\mathit{E}}_{\mathrm{roll}}}{{\mathit{\rho}}_{\mathrm{0}}{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}}}\frac{\mathit{\rho}}{{\mathit{\rho}}_{\mathrm{0}}}\mathrm{~}\sqrt{\frac{{\mathit{E}}_{\mathrm{roll}}}{{\mathit{m}}_{\mathrm{0}}}}\mathit{\phi}\mathit{.}& & \end{array}$(27)where we use Eq. (25). Using the rolling energy of ice particles, c_{s,eff} is given by $\begin{array}{ccc}{\mathit{c}}_{\mathrm{s}\mathit{,}\mathrm{eff}}\mathrm{~}\mathrm{1.1}\mathrm{\times}{\mathrm{10}}^{\mathrm{3}}\mathit{\phi}\mathrm{cm}\mathit{/}\mathrm{s}\mathit{.}& & \end{array}$(28)Therefore, in the case of C_{v} = 3 × 10^{7}, v_{d} is not low enough in the beginning of the simulation, where the aggregate has a low filling factor. However, the boundary velocity difference reaches lower than the effective sound speed when φ ≳ 10^{2}.
3.3. Dependence on the size of the initial BCCA cluster
To confirm that Eq. (25) is valid in the lower density region, we perform the simulations with the different number of particles, which is equivalent to the different sizes of the initial dust aggregates. Figure 6 shows dependence on the number of particles of the initial BCCA cluster.
Fig. 6 Pressure P in [Pa] against filling factor φ with a different number of particles N. Each line shows the average of ten runs of the fixed number of particles: N = 1024,4096, and 16 384. The other parameters are C_{v} = 3 × 10^{7}, k_{n} = 0.01, and ξ_{crit} = 8 Å in the case of N = 1024, 4096, and C_{v} = 1 × 10^{7}, k_{n} = 0.01, and ξ_{crit} = 8 Å in the case of N = 16 384. The dashed line is Eq. (25). 
The initial numbers of particles are 1024, 4096, and 16 384. The other parameters are C_{v} = 3 × 10^{7}, k_{n} = 0.01, and ξ_{crit} = 8 Å in the case of N = 1024 and N = 4096, and C_{v} = 1 × 10^{7}, k_{n} = 0.01, and ξ_{crit} = 8 Å in the case of N = 16 384. We chose lower C_{v} in the case of N = 16 384 to investigate the strength in lower φ region. Each line represents the average of ten runs for each simulation as in Figs. 4b and 5. We draw the averaged line from lower φ than in Fig. 5. In such a low φ region, we consider that the pressure is zero for some runs because the aggregate is isolated from the copies of the aggregate over the periodic boundaries. Except for the initial deviation in low φ, all lines show good agreement with Eq. (25) where φ ≲ 0.1. The result agrees in lower φ for runs with larger N. Therefore, we conclude that the formula Eq. (25) is valid for φ ≲ 0.1.
3.4. Dependence on the normal damping force
As described in Sect. 2.2, we adopt the normal damping force to reduce the normal oscillations in addition to Wada et al. (2007). To confirm that this damping factor does not affect the simulation results, we set the damping factor k_{n} as a parameter. Figure 7 shows dependence of pressure on the normal damping factor k_{n}.
Fig. 7 Pressure P in [Pa] against filling factor φ with different normal damping force. We put the same ten initial conditions varying the normal damping force with k_{n} = 0, k_{n} = 10^{2}, and k_{n} = 10^{1}. Each line shows the result of one run. The other parameters are N = 16 384, C_{v} = 3 × 10^{7}, and ξ_{crit} = 8 Å. 
The fixed parameters are N = 16 384 C_{v} = 3 × 10^{7}, and ξ_{crit} = 8 Å. Each line represents the result of one run for k_{n} = 0, 10^{2}, and 10^{1}, respectively. This figure clearly shows that the normal damping force does not affect the simulation results.
As mentioned in Sect. 3.1, the compressive strength in the lowdensity region (φ ≲ 0.1) is expected to be determined by the rolling motion. To confirm this, we calculate the total energy dissipations of all motions, which are normal damping, rolling, sliding, and twisting. Figure 8 shows the dissipated energy for each mechanism, the dissipated energies in the case without the normal damping and those in the case of k_{n} = 0.01.
Fig. 8 Energy dissipation of each dissipation mechanism in [erg] against filling factor φ. The solid lines show the result in the case without the normal damping and the dashed lines in the case of k_{n} = 0.01 and The results in the case of k_{n} = 10 are not plotted because they are the same as those in the case of k_{n} = 0.01 and indistinguishable. The dissipation mechanisms are normal damping, rolling, sliding, and twisting. The dissipation energy by sliding motion is less than 10^{9} erg. 
The dissipated energy in the case of k_{n} = 10 is indistinguishable from those in the case of k_{n} = 0.01, and thus we do not plot them. The dissipation energy of the sliding force is less than 10^{9} erg, so it is not depicted in this figure. The dissipation by the rolling and twisting is almost the same in the cases with and without the normal damping. Thus, we confirm that the normal damping does not affect the compressive strength, although it dissipates the energy of the normal oscillations. Aside from the normal dissipation, the dominant dissipation mechanism is the rolling motion. This clearly shows that the static compression is determined by the rolling motion of each connection, as mentioned in Sect. 3.1. Where φ ≳ 0.1, the energy dissipation by twisting motion occurs. This is why Eq. (25) is valid until the filling factor reaches 0.1 as mentioned in Sect. 3.1. In the highdensity region, where φ ≳ 0.1, another formulation is required but that is beyond the scope of this paper.
3.5. Dependence on the rolling energy
We also investigate the dependence of the compressive strength on ξ_{crit}. Since E_{roll} is proportional to ξ_{crit}, we investigate the dependence on the rolling energy in this section. Figure 9 shows the dependency on ξ_{crit}.
Fig. 9 Pressure P in [Pa] against filling factor φ with different critical displacement, ξ_{crit}. We put the same ten initial conditions varying ξ_{crit} with ξ_{crit} = 32,16,8,4, and 2 Å, respectively. Each line shows the average of ten runs. The other parameters are N = 16 384 C_{v} = 3 × 10^{7}, and k_{n} = 10^{2}. 
We vary ξ_{crit} with 32, 16, 8, 4, and 2 Å. The fixed parameters are N = 16 384, C_{v} = 3 × 10^{7}, and k_{n} = 10^{2}. This result shows that the compressive strength is almost the same in the lowdensity region. This is because the periodic boundary creates the additional voids as discussed in Sect. 3.1, so we should not focus on the lowdensity region. The lines in the case of ξ_{crit} = 2,4, and 8 Å are on the corresponding lines of Eq. (25) where φ ≲ 0.1. The line in the case of ξ_{crit} = 16 Å has a little deviation, and in the case of ξ_{crit} = 32 Å it has a deviation from their corresponding lines of Eq. (25). The reason the lines in the case of ξ_{crit} = 16,32 Å deviate from the corresponding lines of Eq. (25) is that the dissipation energy is dominated not by rolling motion but by a twisting motion as indicated in Fig. 10. This figure shows that dissipated energy of each dissipation mechanism.
Fig. 10 Energy dissipation of each dissipation mechanism in [erg] against filling factor φ. Each panel represents the case of different ξ_{crit}, which are 32, 16, and 8 Å, corresponding to Fig. 9. 
We show the results of the cases with ξ_{crit} = 8, 16, and 32 Å. The normal damping does not contribute to the compressive strength as discussed in Sect. 3.4, and thus we focus on the rolling and twisting motions.
When ξ_{crit} ≤ 8 Å, the dissipation energy is dominated by rolling motion. For ξ_{crit} = 32 Å, on the other hand, the dissipation energy is dominated by a twisting motion. For ξ_{crit} = 16 Å, the dissipation energy of rolling and twisting motion is comparable, and making it the marginal case. Thus, the reason Eq. (25) is not valid when ξ_{crit} ≥ 16 Å is that the twisting motion is the dominant mechanism for determining the compressive strength. Therefore, we conclude that Eq. (25) is valid when ξ_{crit} ≤ 8 Å.
3.6. Fractal structure
We also investigate how the fractal structure of the dust aggregate changes. Figure 11 shows how many particles are inside the distance r_{in} for four snapshots.
Fig. 11 Number of particles inside the radius r_{in} against normalized radius r_{in}/r_{0}. This figure represents the fractal structure of the compressed aggregates in our simulation for various φ. We set a particle as an origin and count the number of particles inside r < r_{in}, where r is the distance from the origin to each particle’s center. Then we count the same correlation of all particles as an origin and take their average (similar figure of Fig. 7 in the paper of Wada et al. 2008). Each line shows the result at the different time steps. The solid thick lines represents the structure of fractal dimension D = 2, and dashed lines represent D = 3 for each corresponding φ. The dotted line shows the number of particles in calculation. The region below this line corresponds to inside the periodic boundaries. 
We select one run from the case with N = 16 384, C_{v} = 3 × 10^{7}, k_{n} = 10^{2}, and ξ_{crit} = 8 Å. Each snapshot is when φ = 0.003,0.01,0.03 and 0.1, respectively. We take a particle as an origin and count the number of particles inside r < r_{in}, where r is the length from the origin. Then we set the computational region as an origin for all the other particles inside and take an average of them. We obtain the same trend in several runs in the cases of different shapes of initial aggregates.
We also count particles beyond the periodic boundaries. In high r_{in}, $\mathit{N}\mathrm{\propto}{\mathit{r}}_{\mathrm{in}}^{\mathrm{3}}$ because of copies over the periodic boundary distributed as a fractal dimension of three. Therefore, where N(r < r_{in}) ≳ 16 384, N must be $\mathit{N}\mathrm{\propto}{\mathit{r}}_{\mathrm{in}}^{\mathrm{3}}$. However, it is almost out of the range of Fig. 11. Figure 11 shows the number of particles in calculation, which is N = 16 384. The results over this line are affected by the periodic boundary condition and those below this line are in computational region. Thus, the results below the line represent the fractal structure inside the computational region and are not the artificial effect of the periodic boundary condition.
Since the initial aggregate is a BCCA cluster, N is proportional to ${\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}$. In the case of φ = 0.003, which is equivalent to φ of the initial BCCA cluster, $\mathit{N}\mathrm{\propto}{\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}$ as shown in Fig. 11. When the fractal dimension is three, N can be written as $\begin{array}{ccc}\mathit{N}\mathrm{(}\mathit{r}\mathit{<}{\mathit{r}}_{\mathrm{in}}\mathrm{)}\mathrm{=}\frac{\mathit{\phi V}\mathrm{(}\mathit{r}\mathit{<}{\mathit{r}}_{\mathrm{in}}\mathrm{)}}{{\mathit{V}}_{\mathrm{0}}}\mathrm{=}\mathit{\phi}{\left(\frac{{\mathit{r}}_{\mathrm{in}}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathrm{3}}\mathit{,}& & \end{array}$(29)where $\mathit{V}\mathrm{(}\mathit{r}\mathit{<}{\mathit{r}}_{\mathrm{in}}\mathrm{)}\mathrm{=}\mathrm{(}\mathrm{4}\mathit{/}\mathrm{3}\mathrm{)}\mathit{\pi}{\mathit{r}}_{\mathrm{in}}^{\mathrm{3}}$. We also plot this equation for each φ in Fig. 11, with good agreement on a large scale, while maintaining $\mathit{N}\mathrm{\propto}{\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}$ on a small scale.
Therefore, the structure evolution in the static compression is as follows. Initially, $\mathit{N}\mathrm{\propto}{\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}$ because the aggregate is a BCCA cluster. As compression proceeds, the fractal dimension D becomes three on a large scale, while it is two on a small scale. The transit scale from D = 2 to D = 3 becomes smaller as compression proceeds until D = 3 on any scale. This structure evolution means that the static compression reconstructs the aggregate first on a large scale when keeping the smallscale BCCA structure. This is the reason the rolling motion determines the compressive strength, as discussed in Sect. 4.
3.7. Silicate case: comparison with previous studies
The compressive strength has been investigated in the previous study (Seizinger et al. 2012). To investigate the connection of compressive strength from the lowdensity to the highdensity regions, we perform simulations in the case of silicate with the same parameters of Seizinger et al. (2012). Figure 12 shows compression in the case of silicate whose monomer size is 0.6 μm.
Fig. 12 Pressure P in [Pa] against filling factor φ. This figure is same as Fig. 4 but for the case of silicate particles (r_{0} = 0.6 μm). 
The parameters are N = 16 384, C_{v} = 3 × 10^{7}, and k_{n} = 0.01. Figure 12a shows the results of ten runs with different initial aggregates and Fig. 12b shows their average. Using the rolling energy of silicate, which is E_{roll} = 1.42 × 10^{8} erg, we also plot the line of Eq. (25) in Fig. 12b. Since t_{0} is given by 1.71 × 10^{9} s for silicate aggregates, v_{d} becomes 4.01 cm/s for φ = 10^{2} with C_{v} = 3 × 10^{7}. This v_{d} is larger than c_{s,eff}( = 0.77 cm/s when φ = 0.01) for silicate aggregates, allowing the numerical results shown in Fig. 12 to deviate from the line of Eq. (25) in the low φ region. When v_{d} = c_{s,eff}, φ = 3.4 × 10^{2}, the compressive strength should obey Eq. (25) when φ ≳ 3.4 × 10^{2}. In the case of silicate, computational time is huge compared with ice particles. We take a relatively high value of the boundary speed to save on computational time. Therefore, the result deviates from Eq. (25) in the lowdensity region because of the high velocity. In other words, the compression is not static in the low density region. In the highdensity region, on the other hand, the result is in good agreement with Eq. (25), suggesting that Eq. (25) is applicable to aggregates consisting of silicate particles with different r_{0}.
Figure 13 compares our simulation results and Eq. (25) in the lowdensity region (φ < 0.1) with the results of Seizinger et al. (2012) and the fitting formula to experiments (Güttler et al. 2009). This figure corresponds to Fig. 4 in Seizinger et al. (2012). They performed similar Nbody simulations to ours but using a BPCA aggregate composed of silicate particles as an initial condition. The compressive strength of our simulations shows good agreement with the same interaction model in Seizinger et al. (2012) with a little discrepancy: φ = 0.24 at P = 300 Pa in our simulations and φ = 0.21 at P = 300 Pa in Seizinger et al. (2012). The discrepancy, 13% in φ, may be caused by the difference in the initial aggregate or the pressure measurement method. The fitting formula of Güttler et al. (2009) suggests φ = 0.17 at P = 300 in the experiments. The discrepancy from our simulations is 29% in φ. In applicable uses of the static compression formula, we focus on obtaining φ with a given P.
4. Understanding the compressive strength formula
In this section, we analytically derive the compressive strength and confirm Eq. (25). First, we consider the structure of a fluffy aggregate in static compression in our simulations. As described in Sect. 2.3, we adopt the periodic boundary condition and put a BCCA cluster as the initial condition. This corresponds to a large aggregate that filled up with BCCA clusters three dimensionally. As compression proceeds, the initial BCCA cluster is compressed but the aggregate keeps smaller BCCA structure as confirmed in Sect. 3.6. Therefore, the aggregate in static compression always consists of BCCA clusters on some scale and filled up with them. Figure 14 illustrates the aggregate in static compression. The enclosed lines depict BCCA clusters on a small scale.
Fig. 13 The filling factor φ against pressure P in [Pa]. This figure is same as Fig. 12, but plotted with a linear scale of φ and reversal of xy axis to compare with previous studies (see Fig. 4 in Seizinger et al. 2012). The dotted line is the result of numerical simulations in the highdensity region (φ ≳ 0.1) in Seizinger et al. (2012) and the thin solid line is the fitting formula proposed by Güttler et al. (2009). Our results consistently connect to the previous simulations in the highdensity region. 
Fig. 14 Schematic drawing of compression of a dust aggregate consisting of a number of BCCA clusters. The left figure shows a dust aggregate consisting of many BCCA clusters and the BCCA clusters are distributed three dimensionally. Each enclosed line represents each region dominated by the BCCA clusters. The central figure is a BCCA cluster, receiving pressure from the next clusters. The BCCA cluster has a large void depicted in the central figure, and thus the void would be compressed, as expressed in the right figure. The required energy to compress the void is the energy to rotate the connection of monomers in contact. Therefore, the compression can be determined by the rolling motion of monomer connection on the connecting point of the subclusters. 
Next, we consider why the compressive strength can be determined by the rolling energy. The internal mass density and the volume filling factor of the aggregate are equal to those of the BCCA clusters. Compression of the whole aggregate proceeds by compression of each cluster. Therefore, the compressive strength of the whole aggregate would be determined by BCCA clusters. The righthand panel of Fig. 14 illustrates compression of one of BCCA clusters. The pressure on the BCCA cluster is exerted by neighbor clusters, which causes the compression of the BCCA cluster. The BCCA cluster can be divided further into two smaller subclusters because BCCA clusters are created by clustercluster aggregation. A large void exists between the two smaller clusters, and they are connected with one connection of monomers in contact. The compression of the BCCA cluster occurs by crashing the large void, which only requires rolling the monomers at the connection.
We now estimate the compressive strength. In static compression, the aggregate is compressed by external pressure. Each BCCA cluster feels a similar pressure, P. Using the pressure, the force on the BCCA cluster is approximately given by $\begin{array}{ccc}\mathit{F}\mathrm{~}\mathit{P}\mathrm{\xb7}{\mathit{r}}_{\mathrm{BCCA}}^{\mathrm{2}}\mathit{.}& & \end{array}$(30)Since the crashing of the large void is accompanied by the rolling of a pair of monomers in contact, the work required for the crashing is given by socalled the rolling energy of monomers, E_{roll} (Dominik & Tielens (1997) or see Eq. (1) for its definition). Therefore, the required force to compress the aggregate satisfies $\begin{array}{ccc}\mathit{F}\mathrm{\xb7}{\mathit{r}}_{\mathrm{BCCA}}\mathrm{~}{\mathit{E}}_{\mathrm{roll}}\mathit{.}& & \end{array}$(31)Substituting Eq. (30), we further obtain the required pressure to compress the aggregate as $\begin{array}{ccc}\mathit{P}\mathrm{~}\frac{{\mathit{E}}_{\mathrm{roll}}}{{\mathit{r}}_{\mathrm{BCCA}}^{\mathrm{3}}}\mathrm{\xb7}& & \end{array}$(32)The radius of the BCCA clusters can be written by using the physical values of the whole aggregate. The internal density of the BCCA cluster is dependent on its radius. The BCCA cluster has the fractal dimension of two, and its radius is approximately given by r_{BCCA} = N^{1/2}r_{0}, where N is the number of constituent monomers in the BCCA subcluster. The internal density of the BCCA cluster is evaluated as $\begin{array}{ccc}\mathit{\rho}\mathrm{~}\frac{\mathit{N}{\mathit{m}}_{\mathrm{0}}}{{\mathit{r}}_{\mathrm{BCCA}}^{\mathrm{3}}}\mathrm{~}{\left(\frac{{\mathit{r}}_{\mathrm{BCCA}}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{1}{\mathit{\rho}}_{\mathrm{0}}\mathrm{\xb7}& & \end{array}$(33)Using Eqs. (32) and (33), we finally obtain the required pressure (or the compressive strength) as $\begin{array}{ccc}\mathit{P}\mathrm{~}\frac{{\mathit{E}}_{\mathrm{roll}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}}{\left(\frac{\mathit{\rho}}{{\mathit{\rho}}_{\mathrm{0}}}\right)}^{\mathrm{3}}\mathrm{\xb7}& & \end{array}$(34)This is the same as Eq. (25) obtained from our numerical simulations.
5. Summary
We investigated the static compressive strength of highly porous dust aggregates, whose filling factor φ is lower than 0.1. We performed numerical Nbody simulations of static compression of highly porous dust aggregates. The initial dust aggregate is assumed to be a BCCA cluster. The particleparticle interaction model is based on Dominik & Tielens (1997) and Wada et al. (2007). We introduced a new method for compression and adopted the periodic boundary condition in order to compress the dust aggregate uniformly and naturally. Because of the periodic boundary condition, the dust aggregate in the computational region represents one part of a large aggregate, and thus we could investigate the compression of a large aggregate. The periodic boundaries move toward the center, and the distance between the boundaries becomes small. To measure the pressure of the aggregate, we adopted a similar manner to the one used in moleculardynamics simulations. As a result of the numerical simulations, our main findings are as follows.

The relation between the compressionpressure P and the filling factor φ can be written as $\begin{array}{ccc}\mathit{P}\mathrm{=}\frac{{\mathit{E}}_{\mathrm{roll}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}}{\mathit{\phi}}^{\mathrm{3}}\mathit{,}& & \end{array}$(35)where E_{roll} is the rolling energy of monomer particles and r_{0} the monomer radius. We defined the filling factor as φ = ρ/ρ_{0}, where ρ is the mass density of the whole aggregate, and ρ_{0} is the material mass density. Equation (35) is independent of the numerical parameters: the number of particles, the size of the initial BCCA cluster, the boundary speed, and the normal damping force. We confirmed that Eq. (35) is applicable in different E_{roll} and r_{0}. We also analytically confirmed Eq. (35).

Equation (35) is valid where φ ≲ 0.1 in the highdensity region. In the lowdensity region, we confirmed that Eq. (35) is valid for φ ≳ 10^{3} in the case of N = 16 384. From the results of different initial sizes of the aggregates, Eq. (35) is valid in the lowerdensity region in the case of the larger aggregates.

The initial BCCA cluster has a fractal dimension of two in the radius of the cluster, although the whole aggregate has a fractal dimension of three because of the periodic boundary. As compression proceeds, the fractal dimension inside the radius of the initial BCCA cluster becomes three, while the fractal dimension on a smaller scale keeps being two. This means that the initial set up, which is that the fractal dimension on a large scale is three and that on a small scale is two, reproduce the structure of a dust aggregate well in static compression as a consequence. This also supports the compressive strength being determined by BCCA structure on a small scale as shown in the analytical approach.

The static compression in the highdensity region (φ ≳ 0.1) has been investigated in the silicate case in previous studies (e.g., Seizinger et al. 2012). We also performed the numerical simulations in the silicate case and confirmed that our results are consistent with those of previous studies in the highdensity region.
The compressive strength formula allowed us to study how static compression affects the porosity evolution of dust aggregates in protoplanetary disks. In applications to dust compression in protoplanetary disks, we use the compressive strength formula to obtain φ with a given P. Moreover, the obtained compressive strength would be applicable to SPH simulations of dust collisions (Geretshauser et al. 2010, 2011). This application of the static compression process is important future work. In this work, we did not study shear or tensile strengths, but they are also worth investigating in future work.
Acknowledgments
We thank Kohji Tomisaka and Hiroki Senshu for fruitful discussions. We appreciate the careful reading and fruitful comments by the anonymous referee and the editor, Tristan Guillot. A.K. is supported by the Research Fellowship from the Japan Society for the Promotion of Science (JSPS) for Young Scientists (24·2120). S.O. acknowledges support by GrantsinAid for
JSPS Fellows (22·7006) from MEXT of Japan. K.W. acknowledges support by Grantsin Aid for Scientific Research (24540459) from MEXT of Japan.
References
 A’Hearn, M. F. 2011, ARA&A, 49, 281 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J. 2004, in Astrophysics of Dust, eds. A. N. Witt, G. C. Clayton, & B. T. Draine, ASP Conf. Ser., 369, 309 [Google Scholar]
 Blum, J., & Wurm, G. 2000, Icarus, 143, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647 [NASA ADS] [CrossRef] [Google Scholar]
 Geretshauser, R. J., Speith, R., Güttler, C., Krause, M., & Blum, J. 2010, A&A, 513, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Geretshauser, R. J., Meru, F., Speith, R., & Kley, W. 2011, A&A, 531, A166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Greenwood, J., & Johnson, K. 2006, J. Colloid Interf. Sci., 296, 284 [CrossRef] [Google Scholar]
 Güttler, C., Krause, M., Geretshauser, R. J., Speith, R., & Blum, J. 2009, ApJ, 701, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Haile, J. M. 1992, Molecular Dynamics Simulation: Elementary Methods, 1st edn. (New York, NY, USA: John Wiley & Sons, Inc.) [Google Scholar]
 Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, in Protostars and Planets II, eds. D. C. Black, & M. S. Matthews, 1100 [Google Scholar]
 Heim, L.O., Blum, J., Preuss, M., & Butt, H.J. 1999, Phys. Rev. Lett., 83, 3328 [NASA ADS] [CrossRef] [Google Scholar]
 Kempf, S., Pfalzner, S., & Henning, T. K. 1999, Icarus, 141, 388 [NASA ADS] [CrossRef] [Google Scholar]
 Krause, M., & Blum, J. 2004, Phys. Rev. Lett., 93, 1103 [NASA ADS] [Google Scholar]
 Meakin, P. 1991, Rev. Geophys., 29, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Nakagawa, Y., Nakazawa, K., & Hayashi, C. 1981, Icarus, 45, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., Tanaka, H., & Sakagami, M.A. 2009, ApJ, 707, 1247 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paszun, D., & Dominik, C. 2006, Icarus, 182, 274 [NASA ADS] [CrossRef] [Google Scholar]
 Paszun, D., & Dominik, C. 2008, A&A, 484, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paszun, D., & Dominik, C. 2009, A&A, 507, 1023 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seizinger, A., Speith, R., & Kley, W. 2012, A&A, 541, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Smirnov, B. M. 1990, Phys. Rep., 188, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Suyama, T., Wada, K., & Tanaka, H. 2008, ApJ, 684, 1310 [NASA ADS] [CrossRef] [Google Scholar]
 Suyama, T., Wada, K., Tanaka, H., & Okuzumi, S. 2012, ApJ, 753, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, H., Himeno, Y., & Ida, S. 2005, ApJ, 625, 414 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, H., Wada, K., Suyama, T., & Okuzumi, S. 2012, Prog. Theor. Phys. Suppl., 195, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2007, ApJ, 661, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2008, ApJ, 677, 1296 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J., & Cuzzi, J. N. 1993, in Protostars and Planets III, eds. E. H. Levy, & J. I. Lunine, 1031 [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zsom, A., Ormel, C. W., Dullemond, C. P., & Henning, T. 2011, A&A, 534, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1 Schematic drawing of the periodic boundary condition. Each box illustrates a boundary box with a side length L for all directions. When the boundary starts to get closer, the aggregate sticks to the neighboring aggregates over the boundary and is compressed by them. It should be noted that this picture is illustrated in the 2D direction, but our simulations are performed in 3D. 

In the text 
Fig. 2 Schematic drawing to illustrate how the particle velocity is calculated when a particle crosses a periodic boundary. For simplicity, we consider this situation in a 2D field, but we actually calculate this in a 3D situation. We consider that a dust particle is close to the boundary in the left figure. In the next time step, the particle crosses the boundary (dashed circle in the right figure). We put the particle on the other side of the boundary as expressed in Eqs. (8) and (10). The velocity component is converted as expressed in Eqs. (9) and (11). This treatment reproduces the isotropic compression in the velocity field well. 

In the text 
Fig. 3 Snapshots of the evolution of an aggregate under compression in the case of N = 16 384. The top three figures are 3D visualizations. They have the same scale with different time epochs. The white particles are inside a box enclosed by the periodic boundaries. The yellow particles are in neighboring boxes to the box of white particles. For visualization, we do not draw the copies on the back and front sides of the boundaries but only 8 copies of the white particles across the boundaries. Each bottom figure represents projected positions onto 2D plane of all particles in each corresponding top figure. The gray points in the bottom figures correspond to the positions of the white particles in the top figures, and the yellow points correspond to those of the yellow particles in the top figures. Scales are in μm. 

In the text 
Fig. 4 a) Pressure P in [Pa] against filling factor φ. The ten thin solid lines show the results for the initial BCCA clusters with different initial random numbers and thick solid line shows the arithmetic average of the ten runs. b) Pressure P in [Pa] against filling factor φ. Same as the thick solid line in a) plotted with a dotted line of Eq. (25). The parameters are N = 16 384, C_{v} = 3 × 10^{7}, k_{n} = 0.01, and ξ_{crit} = 8 Å. 

In the text 
Fig. 5 Pressure P in [Pa] against filling factor φ with different strain rate parameter C_{v}. Each line shows the average of ten runs of the fixed strain rate: C_{v} = 1 × 10^{7}, 3 × 10^{7}, 1 × 10^{6}, 3 × 10^{6}, 1 × 10^{5}. The other parameters are the same for every ten runs: N = 16 384, k_{n} = 0.01, and ξ_{crit} = 8 Å. The dashed line is Eq. (25). 

In the text 
Fig. 6 Pressure P in [Pa] against filling factor φ with a different number of particles N. Each line shows the average of ten runs of the fixed number of particles: N = 1024,4096, and 16 384. The other parameters are C_{v} = 3 × 10^{7}, k_{n} = 0.01, and ξ_{crit} = 8 Å in the case of N = 1024, 4096, and C_{v} = 1 × 10^{7}, k_{n} = 0.01, and ξ_{crit} = 8 Å in the case of N = 16 384. The dashed line is Eq. (25). 

In the text 
Fig. 7 Pressure P in [Pa] against filling factor φ with different normal damping force. We put the same ten initial conditions varying the normal damping force with k_{n} = 0, k_{n} = 10^{2}, and k_{n} = 10^{1}. Each line shows the result of one run. The other parameters are N = 16 384, C_{v} = 3 × 10^{7}, and ξ_{crit} = 8 Å. 

In the text 
Fig. 8 Energy dissipation of each dissipation mechanism in [erg] against filling factor φ. The solid lines show the result in the case without the normal damping and the dashed lines in the case of k_{n} = 0.01 and The results in the case of k_{n} = 10 are not plotted because they are the same as those in the case of k_{n} = 0.01 and indistinguishable. The dissipation mechanisms are normal damping, rolling, sliding, and twisting. The dissipation energy by sliding motion is less than 10^{9} erg. 

In the text 
Fig. 9 Pressure P in [Pa] against filling factor φ with different critical displacement, ξ_{crit}. We put the same ten initial conditions varying ξ_{crit} with ξ_{crit} = 32,16,8,4, and 2 Å, respectively. Each line shows the average of ten runs. The other parameters are N = 16 384 C_{v} = 3 × 10^{7}, and k_{n} = 10^{2}. 

In the text 
Fig. 10 Energy dissipation of each dissipation mechanism in [erg] against filling factor φ. Each panel represents the case of different ξ_{crit}, which are 32, 16, and 8 Å, corresponding to Fig. 9. 

In the text 
Fig. 11 Number of particles inside the radius r_{in} against normalized radius r_{in}/r_{0}. This figure represents the fractal structure of the compressed aggregates in our simulation for various φ. We set a particle as an origin and count the number of particles inside r < r_{in}, where r is the distance from the origin to each particle’s center. Then we count the same correlation of all particles as an origin and take their average (similar figure of Fig. 7 in the paper of Wada et al. 2008). Each line shows the result at the different time steps. The solid thick lines represents the structure of fractal dimension D = 2, and dashed lines represent D = 3 for each corresponding φ. The dotted line shows the number of particles in calculation. The region below this line corresponds to inside the periodic boundaries. 

In the text 
Fig. 12 Pressure P in [Pa] against filling factor φ. This figure is same as Fig. 4 but for the case of silicate particles (r_{0} = 0.6 μm). 

In the text 
Fig. 13 The filling factor φ against pressure P in [Pa]. This figure is same as Fig. 12, but plotted with a linear scale of φ and reversal of xy axis to compare with previous studies (see Fig. 4 in Seizinger et al. 2012). The dotted line is the result of numerical simulations in the highdensity region (φ ≳ 0.1) in Seizinger et al. (2012) and the thin solid line is the fitting formula proposed by Güttler et al. (2009). Our results consistently connect to the previous simulations in the highdensity region. 

In the text 
Fig. 14 Schematic drawing of compression of a dust aggregate consisting of a number of BCCA clusters. The left figure shows a dust aggregate consisting of many BCCA clusters and the BCCA clusters are distributed three dimensionally. Each enclosed line represents each region dominated by the BCCA clusters. The central figure is a BCCA cluster, receiving pressure from the next clusters. The BCCA cluster has a large void depicted in the central figure, and thus the void would be compressed, as expressed in the right figure. The required energy to compress the void is the energy to rotate the connection of monomers in contact. Therefore, the compression can be determined by the rolling motion of monomer connection on the connecting point of the subclusters. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.