# An even-load-distribution design for composite bolted joints using a novel
circuit model and neural network

Cheng Qiu  Department of Mechanical and Aerospace Engineering, Hong Kong
University of Science and Technology, Hong Kong SAR, China  Yuzi Han
Department of Mechanical and Aerospace Engineering, Hong Kong University of
Science and Technology, Hong Kong SAR, China  Logesh Shanmugam  Department of
Mechanical and Aerospace Engineering, Hong Kong University of Science and
Technology, Hong Kong SAR, China  Fengyang Jiang  School of Aeronautic Science
and Engineering, Beihang University, Beijing, China  Zhidong Guan  Shanyi Du
Jinglei Yang

###### Abstract

Due to the brittle feature of carbon fiber reinforced plastic laminates,
mechanical multi-joint within these composite components shows uneven load
distribution for each bolt, which weaken the strength advantage of composite
laminates. In order to reduce this defect and achieve the goal of even load
distribution in mechanical joints, we propose a machine learning-based
framework as an optimization method. Since that the friction effect has been
proven to be a significant factor in determining bolt load distribution, our
framework aims at providing optimal parameters including bolt-hole clearances
and tightening torques for a minimum unevenness of bolt load. A novel circuit
model is established to generate data samples for the training of artificial
networks at a relatively low computational cost. A database for all the
possible inputs in the design space is built through the machine learning
model. The optimal dataset of clearances and torques provided by the database
is validated by both the finite element method, circuit model, and an
experimental measurement based on the linear superposition principle, which
shows the effectiveness of this general framework for the optimization
problem. Then, our machine learning model is further compared and worked in
collaboration with commonly used optimization algorithms, which shows a
potential of greatly increasing computational efficiency for the inverse
design problem.

##  1 Introduction

Composite laminates, especially carbon fiber reinforced plastics (CFRP), has
been widely used in aerospace industries, due to the excellent strength-to-
weight ratio and high stiffness-to-weight ratio. Although integral-forming of
composite components may take full advantages of its load-carrying ability,
joints among the CFRP parts are inevitable where different connection
techniques may be used including mechanically fastened joints, adhesively
bonded joints, and hybrid joints[1]. Mechanically fastened joints, which
dominates among these techniques in their usage in aerospace structure, hold
the advantage of easily-removability, non-sensitivity to the environment, and
high load-carrying ability[2, 3].

Although the design process of metal joints has been well established, it has
to be faced with some differences due to the brittle feature of CFRP
laminates. While the plasticity of metal material renders it the load
redistribution ability, the load distribution ratio is almost unchanged for
CFRP mechanical joint before ultimate failure. Generally, the design
methodology or strength validation of the mechanical joint takes the following
steps[4, 5, 6]. First, the stress conditions of the multi-bolt joint should be
extracted from an overall load analysis on global structures. Second, the load
distribution would be determined based on the coordination relationship among
the deformations of each joined part. After acquiring each bolt load, detailed
stress analysis is going to be performed on the area where damage is most
likely to occur[7, 8].

Determining bolt load distribution is one of the research hotspots in
composite joints since it plays a fundamental role in structural analysis[9,
10, 11]. The typical problem is a joint with single-column, multi-row bolts
under external tensile load. Experimental measurement of bolt load is
difficult as a result of the narrow space at the contact area for the
placement of sensors like strain gauges[12]. Finite element method or
analytical method with engineering experience parameters could be used but
faced with drawbacks of computational efficiency and accuracy[10, 13]. The
spring-based model, which transforms different parts in the joint structure
into springs with specific stiffness, has been widely used in the literature
due to its simplicity to implement[14]. In the meantime, a new requirement for
the spring-based model has been put forward with regard to the importance of
considering the effect of friction, the relevant parameters of which include
bolt-hole clearance and tightening torques. Instead of using a constant spring
stiffness, the load-displacement curve shows a four-stage shape considering
the transition of contact force from static frictions to sliding frictions,
making a variable spring stiffness in the model[15, 9, 16]. While researchers
have clarified influencing factors of bolt load distribution, optimization of
an even-load-distribution design has seldom been focused on.

Machine Learning, which discerns hidden patterns from complex dataset has been
proven to be a very efficient tool in many research fields including solid
mechanics[17, 18, 19, 20]. As a branch of artificial intelligence, on one
hand, it can be incorporated as an alternative way of generating mechanical
response of a given structure or compounds which used to obtain by advanced
computational techniques like density functional theory, molecular dynamics,
and finite element method[21, 22, 23]. On the other hand, it is claimed to be
a promising way to accelerate the material design process which involves an
inverse problem of designing new materials or structures with desired
property[24, 25, 26]. Some interesting and meaningful applications includes
designing hierarchical structures with high fracture toughness[27], predicting
residual strength from acoustic emission data[28], solving inverse structure-
process problem in addictive manufacturing[29].

In this work, we aim to provide an even-load-distribution design process for
multi-bolt CFRP composite joint using a machine learning based framework.
Based on the spring model, we put forward a simple equivalent circuit model in
Matlab/Simulink that can easily produce dataset for the training of machine
learning model. The bolt-hole clearance and tightening torques are chosen as
design inputs due to their significant influence on load distribution as well
as can be easily manipulated in practice. A symbol of load distribution
unevenness is defined as output to represent the level of differences in bolt
loads before failure. The artificial neural network is trained as a response
predictor and then used to generate a database for all the possible inputs in
our design space. Finally, the optimal parameters provided by this general
machine learning-based framework is validated for its reducing the load
distribution unevenness in a three-bolt joint.

##  2 Problem description

In this section, the load distribution problem of a typical single-lap multi-
bolted joint is introduced. The corresponding circuit model using
Matlab/Simulink is established in which bolt-hole clearance and frictions
caused by the tightening torque can be taken into consideration. After a brief
introduction of the circuit model, the effectiveness and accuracy of the model
are validated through a comparison with FEM.

###  2.1 Circuit model for predicting bolt-load distribution

For the load-sharing analysis of multi-bolted joint structure as shown in
Fig.1, the spring-mass model is the most commonly used method which holds the
advantage of low computational cost and high efficiency against other methods
like FEM, analytical method, and experimental measurement. With the increasing
number of bolts, obtaining bolt load from the spring-mass model becomes
complicated especially when various bolt-hole clearances and tightening
torques are considered, which may involve non-linear stages. To solve this
problem, in our previous work[15], a novel circuit model is established for
determining the bolt-load distribution and joint stiffness. Then, in the
circuit model, bolt load distribution is reflected by the flow of current.
Here, a brief introduction of this circuit model is given.

Figure 1: (a) A typical structure of three-bolt, single-lap joint (b) Spring-
mass model of multi-bolt joint

The spring-based method, which is a stiffness method essentially, separating
the global deformation under external load into several extensions of
different parts. These parts are assumed to have a simple load-deformation
relationship which can be represented as springs in the model. The stiffness
of each spring is the key parameter, especially for the $K_{bolt}$ shown in
Fig.1 denoting the load-deformation relationship in the bolt area.

$K_{plate}$ in the spring-based model, which is the tensile stiffness of the
jointed plate, can be easily calculated as,

$K_{plate}=\frac{E_{px}A_{p}}{l_{p}}$ (1)

where $E_{px}$ is the laminate modulus along the applied load direction.
$A_{p}$ stands for cross-sectional area of the plate while $l_{p}$ is the bolt
pitch.

It is revealed that with the introduction of bolt-hole clearance, a four-stage
curve of load-displacement behavior occurs for the single-lap joint. In the
beginning, when the applied load is too small to invoke any sliding, the
external load is counteracted by the static frictions from the upper and lower
contact surfaces on plates(which is highlighted by red and blue color in Fig.1
). In the second stage, as the applied load increases, the contact forces
would exceed the maximum static friction, causing the slides at the contact
area. This sliding is proven to take place first at the contact surfaces
between two jointed plates since that ratio of $F_{1}$ to $F_{2}$ is greater
than 1 shown as,

$R=\frac{F_{2}}{F_{1}}=(\frac{64t^{2}}{3d^{2}(1+\mu_{b})}+2.4)\cdot\frac{A_{0}G_{p}}{A_{b}G_{b}}+1$
(2)

Here $d$, $\mu_{b}$, $A_{b}$ and $G_{b}$ are the diameter, Poisson’s ratio,
minimum cross-sectional area and shear modulus of bolt respectively. Whilst,
$A_{0}$ is the contact area between bolt head and plate. $G_{p}$ is the shear
modulus of laminate.

When the increasing external load reaches the summation of maximum static
frictions on both upper and lower contact surfaces, a global sliding happens
to lead to a plateau phase is shown in the load-displacement curve. Once the
bolt rod contacts with the edge of bolt hole on laminates, the load-
displacement curve goes into the fourth stage in which the joint stiffness
keeps a constant value until final failure.

According to the above analysis, joint stiffness would vary in each of the
four stages as a result of the changes in friction condition, that is,

(1) Phase 1 ($F_{2}<f_{static}$)

$\frac{1}{K_{bolt}^{(1)}}=(\frac{12t}{5A_{b}G_{b}}+\frac{8t^{3}}{3E_{b}I_{b}}+\frac{(3-R)t}{4A_{0}G_{p}})\cdot\frac{1}{1+R}+\frac{(4l_{p}-3\Phi)t^{2}}{32E_{px}I_{p}}$
(3)

(2) Phase 2 ($F_{1}<f_{static}=F_{2}$)

$\frac{1}{K_{bolt}^{(2)}}=(\frac{12t}{5A_{b}G_{b}}+\frac{8t^{3}}{3E_{b}I_{b}}+\frac{3t}{4A_{0}G_{p}})+\frac{(4l_{p}-3\Phi)t^{2}}{32E_{px}I_{p}}$
(4)

(3) Phase 3 ($F_{1}=f_{static}=F_{2}$)

$\frac{1}{K_{bolt}^{(3)}}=\infty$ (5)

(4) Phase 4 ($displacement>clearance$)

$\frac{1}{K_{bolt}^{(4)}}=\frac{4t}{3G_{p}A_{p}}+[\frac{4}{tE_{b}}+\frac{2}{t(\sqrt{E_{px}E_{py}})}](1+3\beta)$
(6)

In these equations, $t$, $\Phi$ are geometric parameters can be found in
Fig.1. $E_{b}I_{b}$ is the bending stiffness of bolt rod. While bending moment
in single-lap joints can be reacted by both laminates and bolts, the term
$\beta$ here represents the fraction of the bending moment which is reacted by
the non-uniform contact stresses in laminates. Here $\beta=0.15$ for a typical
protruded-bolt joint. Detailed explanation for all the above equations can be
found in Ref[15].

After figuring out the corresponding stiffness of each spring, the next step
is to acquire bolt load, global deformation, and global joint stiffness from
spring-based model. In order to solve this problem, a novel circuit model is
proposed using Matlab/Simulink. A basic model for the three-bolt joint is
shown in Fig.2. In this circuit model, the applied voltage equals to the
displacement at the loading end. Then the resistance value of each resistor is
set to be reciprocal of the corresponding spring stiffness expressed in
Eq.1-6, making the flow of current equals to bolt load distribution. The
Matlab/simulink provides a simple way to simulate and run the whole process.
Then, results can be outputted to Matlab workspace for the following post-
processing and visualization work.

Figure 2: The circuit model for bolt load sharing analysis of a single-lap,
three-bolt joint

Special attention should be paid to the “$R_{bolt}$” module. The established
module for bolt stiffness is shown in Fig.3. A total of three conditional
switches are incorporated to evoke the three knee points in the four-stage
load-displacement curve. While the transition conditions are presented in the
form of force for the first two knee points in Eq.3 and Eq.4, the form of
displacement which corresponds to the voltage in circuit model is adopted in
this module. In practice, after measuring the voltage value using a sensor
block, its value is extracted and compared to “a”, “b” and “c” to control the
switches. “a”, “b” and “c” denoting for the three displacements of transition
conditions are defined and imported into the circuit model using a Matlab
script. Finally, the bolt load is measured by a current sensor block and
outputted to Matlab workspace by the output modules shown in Fig.2.

Figure 3: Module for bolt stiffness considering the effect of bolt-hole
clearances and tightening torques

The effect of bolt-hole clearances and tightening torques can be reflected by
the values of “a”, “b” and “c”. Among others, the value of “c” here accounts
for bolt-hole clearance. According to the Coulomb’s friction law, the maximum
static friction is proportional to tightening force which relates to
tightening torques and geometries of the bolt. Thus, the displacement-control
knee points can be calculated as,

$a=\frac{f_{static}(1/R+1)}{K_{bolt}^{(1)}}=\frac{vT}{kd}\cdot\frac{(1/R+1)}{K_{bolt}^{(1)}}$
(7)
$b=a+\frac{f_{static}(1-1/R)}{K_{bolt}^{(2)}}=a+\frac{vT}{kd}\cdot\frac{(1-1/R)}{K_{bolt}^{(2)}}$
(8)

Here, $v$ is the friction coefficient between contact surfaces. $T$ is the
applied tightening torque and $k$ is a constant value accounting the relation
between bolt torque and tightening force. $k=0.2$ and $v=0.3$ is adopted in
this paper according to Ref[30].

###  2.2 Validation of circuit model

A three-bolted joint with its geometries shown in Tab.1 is studied in this
paper. The composite layup is $[+45/0/-45/90]_{3S}$ for both laminates, whose
nominal thickness is 3 mm. The lamina properties are $E_{1}=195.00$ GPa,
$E_{2}=8.58$ GPa, $G_{12}=4.57$ GPa. The homogeneous mechanical properties of
laminate are computed using the Classical Laminate Theory and listed in Tab.1
as well as the mechanical parameters of the bolt.

In the circuit model, the resistors marked as “$R_{plate}$” are assigned with
a value reciprocal to Eq.1 while the resistors in the “$R_{bolt}$” module have
a value that corresponds to Eq.3-Eq.6 respectively. It has to be noted that a
fixed-step solver option is adopted in Matlab/Simulink. Key parameters in
solver control include value of the provided Voltage Source $V_{s}$, single
step time $t_{0}$ and total simulation time $t_{s}$. With these parameters,
the total simulated displacement can be expressed as,

$V_{total}=\frac{t_{s}}{t_{0}}\cdot V_{s}$ (9)

In practice, for each step, the Voltage Source provides a 0.005V voltage
increment. With the combination of solver settings for analysis time 60s and
fixed-step time 0.1s, the total applied voltage is 3V, that is 3mm
displacement applied at the loading end.

A three-dimensional numerical model is established in Abaqus software to
provide FEM results as a reference to our circuit model. All of the five parts
contained in the numerical model are 3D solid with reduced-integration, 8-node
linear brick element C3D8R, as illustrated in Fig.4. A radiation-shape mesh
with a bias ratio of 5 is adopted to refine mesh size near each bolt hole.

While the three bolts are assigned with a homogeneous isotropic material
property, each of the laminate is partitioned into 6 layers along its
thickness in which each layer represents $[+45/0/-45/90]$ ply in laminates.
Then, the lamina properties are assigned to each ply along with a local
material orientation according to its ply angle.

A total of 5 contact pairs are included in the model using a surface-to-
surface contact method as shown in the figure. According to Abaqus user’s
manual[31], surfaces with coarse mesh size and greater material stiffness are
more suitable to set as the master surfaces to prevent penetration between
adjacent surfaces. While the default “Hard” contact is used for the normal
contact property, a penalty method with friction coefficient 0.3 is defined to
simulate tangential contact behavior. All the degrees of freedom of both ends
of the assembled structure are coupled to a controlling point to apply
displacement and constraints.

The pre-tightening force is introduced by using an expansion method. As listed
in Tab.1, the bolts are assigned with three thermodynamic parameters
$\alpha_{1},\alpha_{2},\alpha_{3}$ which represent the coefficient of thermal
expansion along local 1, 2, and 3 direction respectively. The $\alpha_{1}$ and
$\alpha_{2}$ are set to be zero to make sure that there is no other contact
force except for the pre-tightening force throughout this temperature drop
process. It is noted that the value of $\alpha_{3}$ and also the temperature
change are chosen and adjusted to make the normal contact force between
laminates and bolt head to be the desired pre-tightening force. A two-step
calculation is conducted in Abaqus, in which the pre-tightening force is
introduced by a temperature drop in Step 1 while the tensile loading is
applied in Step 2. Here in this model, a temperature drop of $60^{o}$ is
defined in this step to generate a normal contact force of 4.3 kN, which
corresponds to a tightening torque of 7 Nm.

Figure 4: 3D finite element model of a three-bolt joint containing 5 contact pairs Table 1: Geometries and mechanical properties of a three-bolt composite joint Geometries |  $l_{p}$(mm) |  $d$ (mm) |  $\Phi$(mm) |  h(mm) |  width (mm) |
---|---|---|---|---|---|---
Value |  60 |  8 |  14 |  8 |  20 |
Laminate |  $E_{px}$(GPa) |  $E_{py}$(GPa) |  $G_{p}$(GPa) |  |  |
Value |  71.60 |  71.60 |  4.57 |  |  |
Bolt |  $E_{b}$ (GPa) |  $G_{b}$ (GPa) |  $\mu_{b}$ |  $\alpha_{1}(/^{o}C)$ |  $\alpha_{2}(/^{o}C)$ |  $\alpha_{3}(/^{o}C)$
Value |  109.78 |  41.27 |  0.33 |  0 |  0 |  0.0001

First, a three-bolt joint with uniform bolt-hole clearances and tightening
torques is studied by both FEM and our circuit model. The bolt-hole clearance
is 0.1 mm which is acquired by slightly reducing the radius of the bolt. And
the tightening torques are equally to be 7 Nm, acquired by adjusting the
expansion parameters of bolt according to the above analysis.

Comparison between results by FEM and circuit model is shown in Fig.5. A
similar trend is shown by the load-displacement curves of both models
illustrated in Fig.5(a). In the beginning, when the applied load is quite
small, the overall stiffness is determined by the deformation of bolt and
laminates under friction. A good agreement is found between circuit model and
FEM which represents the accuracy of Eq.3 in dealing with structural
deformation under friction. As the load increases, sliding firstly comes out
in bolt 1 and bolt 3, due to the fact that bolt load shared by bolt 1 and bolt
3 is greater than the central bolt from the beginning. Then, this changed
stiffness of bolt 1 and bolt 3 would affect the overall bolt load
distribution, causing more time for bolt 2 to enter into the sliding phase.
While the transition on the load-displacement curve caused by the stiffness
change is a piecewise form for the results of circuit model, a rather smooth
curve in shown by the curve of FEM. One of the reasons may be the more
delicate consideration of contact by FEM such as the change of contact area,
the non-uniform contact force, and the slight rotation of contact direction.
Another reason can be the fewer data points for fitting the curve by FEM. As
the load continues to increase, full contact between bolt rod and hole is
achieved. A good agreement can also be found in the curves of three bolt
loads, as shown in Fig.5(b)-(d).

Bolt load distribution ratio predicted by the circuit model is shown in
Fig.5(e). As a result of the symmetry in three-bolt joint, bolt loads shared
by bolt 1 and bolt 3 are exactly the same, which are greater than load shared
by bolt 2. Several platform stages are undertaken in the beginning due to the
transition of the contact status at the bolt area. After that, it is shown
that in this condition of uniform bolt-hole clearance and tightening torque,
the maximum ratio is 37.2% while the minimum ratio is 25.6%. In this case,
when the total appiled load is 50 kN, a difference of 5.8 kN can take place
between bolts.

Figure 5: Comparison between results by FEM and circuit model (a)Load-
displacement curve of the whole structure (b)(c)(d)Bolt load shared by each
bolt vs. global displacement (e)Load distribution ratio of three bolts vs.
global displacement

###  2.3 Effect of the randomness of clearances and tightening torques

In the last section, the effectiveness of our circuit model for the load
distribution calculation is validated by FEM. Here, the effect of the
randomness of clearances and tightening torques is investigated. Since that
the circuit model is built in Matlab/Simulink, random parameters are produced
by a Matlab script and exported into the circuit model. For the clearances, a
random number ranges from 0 to 2 mm is chosen for each bolt, while the random
tightening torques range from a hand-tight 0.5 Nm to a highly-tighten torque
15 Nm.

In the case shown in Fig.6, the parameters for bolt-hole clearances are
(bhc1=1.04mm, bhc2=0.63mm, bhc3=0.12mm) and the torques are (T1=8.41Nm,
T2=9.25Nm, T3=5.90Nm). Compared to Fig.5 where the parameters are identical,
the bolt load-displacement during the loading process has significantly been
affected. In Fig.5, the load shared by bolt 1 and bolt 3 is equal. With the
randomness of the clearances and torques, this equality is changed. In the
case of Fig.6, bolt-hole clearance bolt 3 is smaller compared to other two
bolts, resulting in a very quick jump into the full contact phase. As stated
before, the changing stiffness of bolt 3 affects the overall load distribution
of joint, causing a different transition point and duration of the plateau on
load-displacement curves of the other two bolts. Moreover, it can be seen that
the total applied force illustrated in Fig.6(a) does not show an obvious
piece-wise shape caused by the friction effect, while in fact a highly uneven
load distribution is attained within joint. Load distribution ratios presented
in Fig.6 show a complicated pattern where the maximum ratio can be up to
64.1%. However, since a constant stiffness is achieved if all three bolts turn
into full contact, a converged value of load ratio could be obtained with
increasing external load. Results suggest that at the displacement of 3mm, the
maximum ratio changes from the beginning value 64.1% to a value around 50%,
while the minimum ratio changes from the beginning value 13% to a final value
around 22%.

Comparison of the bolt load distribution ratio and bolt load (at total load
=30 kN) in the three-bolt joint with identical and random bolt parameters is
shown in Tab.2. It can be seen that, with the varied clearances and torques,
load distribution among these three bolts can be affected. Therefore, this
effect of various parameters inspires us to find a set of parameters that can
help to narrow down the differences in the load distribution.

Table 2: Comparison of bolt load distribution ratio between identical and random bolt parameters under a total load of 30 kN Cases |  Parameters (bhc1, bhc2, bhc3, T1, T2, T3) |  Bolt 1 |  Bolt 2 |  Bolt 3
---|---|---|---|---
Identical |  (0.1, 0.1, 0.1, 7.0, 7.0, 7.0) |  11.2KN (37.2%) |  7.7 kN (25.6%) |  11.2 kN (37.2%)
Random |  (1.04, 0.63, 0.12, 8.41, 9.25, 5.90) |  4.98 kN (16.6%) |  7.94 kN (26.5%) |  17.08 kN (56.9%)
Figure 6: Results of the three-bolt joint with random bolt parameters of
clearances (1.54mm, 0.63mm, 0.12mm) and tightening torques (8.41Nm, 9.25Nm,
14.90Nm) (a)Bolt load-displacement curve (b)Bolt load ratio

##  3 Neural network for even-load-distribution design

In this section, by adopting the back-propagation neural network, we are
trying to form a tool for designing a multi-bolt joint having an even-load-
distribution under the common tensile load.

The general framework of this machine learning-based optimization work can be
summarized by Fig.7. It can be seen that our work can be divided into three
parts, that is “Data Generation”, “Neural Network”, and “Find Data”. First,
the input variables include three bolt-hole clearances and three tightening
torques. As explained before, the design space for bolt-hole clearance is
between 0 to 2 mm and for tightening torque is between 0.5 to 15 Nm. Here, an
increment size of 0.2 mm is used for clearances which makes 1000 possible
combinations for three-bolt clearances. In the meanwhile, an increment size of
0.5 Nm is adopted for torques which makes 27000 possible combinations for
three-bolt torques. Thus, a total of $2.7\times 10^{7}$ possible test patterns
are assumed. The following symbol for the evaluation of unevenness is adopted
as outputs:

$u=\frac{max\\{F_{i}\\}-min\\{F_{i}\\}}{max\\{F_{i}\\}+min\\{F_{i}\\}}\quad(i\
for\ total\ load=30kN)$ (10)

Here, $max\\{F_{i}\\}$ indicates the maximum bolt load when to total applied
load is 30kN while $min\\{F_{i}\\}$ indicates the minimum bolt load. If the
unevenness of load distribution is low, this $u$ value should be small.

In Data Generation, 1000 random inputs and corresponding outputs are generated
by Matlab script using our circuit model. In practice, this step takes about
20 minutes mainly due to the call and calculation of the circuit model in
Matlab/Simulink. Then, these 1000 samples are used for the training process of
our back-propagation neural network. It is noted that the 1000 samples are
divided into training set, validation set, and test set using a ratio of
70%/10%/20%. The architecture of the established neural network is shown in
Fig.7. Four hidden layers are included and the number of neurons in each layer
is 30/40/40/30 respectively. The stochastic gradient descent back-propagation
optimizer “Adam” is used as the learning algorithm while the rectified linear
function (ReLU) is used as the activate function for calculation between
layers. After the training and evaluation of neural network, it is used as a
response predictor to obtain a symbol of the unevenness of the test pattern. A
database is thus formed by these test patterns and their predictions. With
this database, a script of pattern generator is written to get the optimal
parameters of bolt-hole clearances and tightening torques which make the
smallest symbol of distribution unevenness. Finally, the optimal results are
tested and evaluated by comparing with FEM and results by the circuit model.

Figure 7: Machine-learning based framework for a even-load-distribution design
of three-bolt joint (a)Data Generation: using circuit model in Matlab to
generate 1000 samples (b)Neural Network: training and validation of the
established neural network (c)Find Data: using neural network to form a
database and find desired data of minimum symbol of uneven load distribution

The ultimate goal of a neural network is to find a group of weights and bias
for every neurons to make a minimum loss function, which indicates the
difference between predictions and target values. Here, a weighted mean square
error (MSE) is defined as loss function to increase the accuracy of predicting
a relatively low value of unevenness $u$, which is important for our
optimization problem. The adopted loss function is shown in Eq.11. Here,
$y_{i}$ is the predicted value and $\hat{y_{i}}$ is the target value. Compared
with the commonly used MSE, a weighting factor $1/(\hat{y_{i}}+0.001)$ is
added, which will make a larger value of loss function when the target value
$\hat{y_{i}}$ is smaller. With the incorporation of this customized loss
function, the accuracy of our neural network is increased with $\hat{y_{i}}$
in the range (0, 0.1) as presented in Fig.8(a) and (b).

$\text{Weighted
MSE}=\frac{1}{n}\sum_{i=1}^{n}\frac{(y_{i}-\hat{y_{i}})^{2}}{\hat{y_{i}}+0.001}$
(11)

Having built the neural network architecture using a widely used Python
machine-learning package “TensorFlow”, we then focus on the validation and
accuracy of the response predictor. One common way to evaluate the performance
of this regression model is through the loss function, i.e. the weighted MSE.
A smaller loss function value means better performance. Fig.8(c) shows the
history of the weighted MSE with increasing epoch number. Every time the epoch
number increases, the whole dataset is used and the regression parameters
within neural network are updated according to the chosen learning algorithm,
thus making a smaller MSE if they are updated towards the right direction. It
can be seen that with the increasing epoch number, MSE gradually decreases
from 0.5 to the minimum MSE 0.0015. In order to avoid over-fitting, we set a
stopping criterion that training process ends when the loss function MSE does
not improve more than 1e-5 for 20 consecutive epochs. So the training process
ends at epoch 355. Another common representation of the performance for this
machine learning model is through a comparison between every single prediction
and its reference value. This comparison is shown in Fig.8(d) in which x-axis
is the reference value and y-axis is the predicted output. A closer position
of a dot to the bisection line ($y=x$) suggests a more accurate predicted
value. In the figure, this is also reflected by the regression coefficient
”$R^{2}$” ranging from 0 to 1. The closer the value to 1, the better matching
of the neural network. It can be seen that the value of R is very close to 1
for both training data and test data. And the band width of the dots is
narrow, indicating a satisfying overall performance of the neural network.

Figure 8: Predicting accuracy with target values in the range (0, 0.1)
(a)normal MSE (b)Weighted MSE; Performance of the neural network model
(c)History of MSE with increasing epoch number (d)Regression coefficient for
training and test dataset

After successfully validated the robustness of neural network, it is used to
generate all the corresponding predictions of all possible $2.7\times 10^{7}$
test patterns. It is noted that the precision of this database is controlled
by the increment size of clearances and torques. Considering the actual
assembly equipment and measurement difficulty, we think that the step size of
0.2 for clearances is acceptable. Since that the effect of tightening torques
is directly reflected by the friction which is nearly one-fifth of the torque
in our cases according to Eq.7, a larger step size of 0.5 for torques is
acceptable. According to this database, the optimal parameters of the desired
property, which is the minimum symbol of unevenness, are as follows:

  * •

Clearances: $bhc1=0.4mm,bhc2=0.2mm,bhc3=0.4mm$

  * •

Torques: $T1=11Nm,T2=10Nm,T3=15Nm$

Since that multi-bolt joint using identical parameters illustrate a relatively
low load ratio shared by bolt 2, the optimal parameters provided by neural
network indicate a smaller bolt-hole clearance at bolt 2 will help to
alleviate the load carried by bolt 1 and 3 while increase the level bolt 2 is
occupied. The detailed reason behind this can be explained through a careful
examination of the loading response of the three-bolt joint with the optimal
parameters, which is shown in Fig.9. It can be seen that joint with these
chosen parameters renders bolt 2 a high initial bolt load ratio, which is up
to 41% when the total applied load is 18.0 kN and displacement is 1.25 mm.
Since that our algorithm only cares about the load distribution when the total
applied load is 30kN, this high unevenness at the beginning phase is not taken
into account. According to the results shown in Fig.5, when all three bolts
enter into the phase of full-contact, the load distributed on bolt 2 is
clearly less than the other two bolts. Therefore, for results of Fig.9, the
high initial load ratio of bolt 2 start to decrease as full contact begins.
Then a low value of unevenness is successfully acquired at the total load of
30kN. The load ratio shared by each bolt is 32.1%/34.7%/33.2% relatively,
which shows a good even-load-distribution.

Figure 9: Loading response of three-bolt joint with optimal parameters for
even-load-distribution design provided by machine learning based framework
(a)Bolt load-displacement curve (b)Bolt load ratio

Although this machine learning-based framework can be successfully implemented
to give a satisfying set of optimal parameters than can fulfill the demand of
minimum unevenness, its efficiency needs to be valued. It takes almost 20
minutes for this whole design campaign of a specific multi-bolt joint. This is
a huge advantages of our circuit model built in Matlab/Simulink over FEM since
our FE model would take about 5 hours generating one single data. The time-
consuming operation of FEM usually comes from the introduction of contact and
three-dimensional elements, which may increase the difficulty of convergence.
A further effort can be made by using parallel computing in Matlab, or decode
the circuit model and replace it by a equal script.

##  4 Discussion

###  4.1 Comparison with other optimization algorithm

In the above work, we use this machine learning model to efficiently build a
database for all the possible combinations of design variables in the given
design space, and then can directly find the best outcomes. A common way to
solve the optimization problem is by using optimization algorithm like Genetic
Algorithm(GA), Ant Colony Optimization(ACO), and Particle Swarm Optimization
(PSO). These advanced algorithms are inspired by the nature for finding the
right iteration direction given initial dataset. Here, results and
computational efficiency of these algorithms, machine learning model, and also
their combinations are compared.

The work flow of the optimization process is shown in Fig.10. Here, both
Genetic Algorithm and Particle Swarm Optimization algorithm provide a way to
approach minimum fitness function value at each iteration. When the stopping
criterion is reached, which is typically evaluating whether the fitness
function has converged, the optimization process stops and best output is
reported. With the circuit model and machine learning model as the predictor
of fitness value, a total of four methods are compared to our machine
learning-based framework, that is, (a) GA+circuit model (b) PSO+circuit model
(c) GA+machine learning model (d) PSO+machine learning model.

For the operation of these algorithms, stopping criterion is one of the most
important factors affecting the computational efficiency. A same stopping
criterion is adopted for all the optimization process, which will end the
progress when the relative change in best objective function value over the
last 20 iterations is less than 1e-3. The optimization history curves of GA
and PSO are shown in Fig.11. It can be seen that both PSO and GA can quickly
reach a convergence at around 13 iterations. The computational time cost
listed in Tab.3 clearly shows that, on one hand, using PSO can have slightly
greater time-efficiency than GA. On other hand, it can be shown that by using
machine learning model as a surrogate model for circuit model, the time cost
by the optimization process can be greatly reduced. For each configuration,
only 1.2 second is needed. This significant improvement is partly because of
the fact that neural network algorithm is quite simple and easy to implement,
and partly because the circuit model-based optimization process needs to call
the Matlab/Simulink module at each iteration. It has to be noted that 20
minutes have already been spent on generating training dataset for this neural
network. However, firstly the training data needed can be further reduced and
secondly, the machine learning model can be expanded when designing joints
with other geometric configurations.

Comparing the optimized parameters provided by these methods and also by our
database using the framework shown in Fig.7, it can be found that the
unevenness values shown in Tab.3 are smaller than the previous results shown
in Fig.9, indicating a better design. This is due to the prescribed precision
in our database as bolt-hole clearance increases by 0.2mm while torque
increases by 0.5 Nm. But using these optimization algorithms, the input
precision can be set to 0.01, which greatly expand the whole design space.
Therefore, better results can be found. Here, the optimized parameters listed
in Tab.3 are different by different methods. However, all of their
corresponding bolt load -displacement curves as presented in Fig.12 show good
evenness of bolt-load-distribution. This suggests that our optimization
problem here may have local minimums, which makes many different design
parameters can lead to a close value of unevenness.

Figure 10: Flowchart of genetic algorithm (GA), particle swarm optimization algorithm (PSO), and the combination with machine learning model Figure 11: The optimization history curves of (a)GA and (b)PSO at each iteration Table 3: Comparison of results and computational efficiency of the optimization methods Methods |  Best input (bhc1, bhc2, bhc3, T1, T2, T3) |  Best output (unevenness u) |  Time
---|---|---|---
GA+Circuit |  (1.0, 0.9, 1.0, 3.2, 10, 1.2) |  0.019 |  86.7min
PSO+Circuit |  (0.32, 0.21, 0.28,6.6, 9.6, 2.2) |  0.007 |  67.8min
GA+ML |  (0.58, 0.41, 0.60, 3.9, 4.2, 12.1) |  0.031 |  9.1s
PSO+ML |  (0.2, 0.2, 0.2, 1.2, 12.6, 0.5) |  0.021 |  1.2s
Figure 12: The bolt load -displacement curves calculated by circuit model
using the optimized parameters provided by the four optimization methods:
(a)GA+Circuit model, (b)PSO+Circuit model, (c)GA+Machine learning, and
(d)PSO+Machine learning

###  4.2 Experimental study for validation

Followed by the intuition given by the optimization results, several composite
three-bolt joints were fabricated for validating the bolt load distribution.
Carbon fiber/epoxy prepregs were used for fabricating, which were cured under
80 degree for 2hours. And then the cured plate was cut by a waterjet cutting
machine ProtoMax with high positional accuracy of 0.05mm. For all the
composite laminates, a layup of $[+45/0/-45/90]_{3S}$ is used. Their
homogeneous laminate property and geometries are listed in Tab.1. The bolt
hole diameters are different for the two groups of test specimens to generate
different bolt-hole clearances when assembled with same steel bolts.

The key problems for this experimental validation is how to evaluate the bolt
load distribution in this multi-bolt structure. Here, a method of calculating
the bolt load ratio using four strain gauges is introduced. The positions of
these four strain gauges on the composite laminates are shown in Fig.13. To
deeply expatiate on the feasibility of the method, a detailed stress analysis
on the bolt hole is performed. For the external load applied to a composite
joint, it can be seen as a combination of bearing load and by-pass load.
Bearing load $P_{br}$ is the load transferred by the bolt while the by-pass
load $P_{by}$ is the remaining load taken by the composite laminates. Due to
the brittle nature of composites, a superposition principle can be used under
linear assumption. According to Ref[32], two non-dimensional stress
concentration factors $\alpha_{br}$ and $\alpha_{by}$, can be utilized to
describe the linear relation between the stress at a specific point and the
applied external stress. Then, the stress at point 1 and point 2, can be
written as,

$\sigma_{1}=\alpha^{1}_{br}\cdot\sigma_{br}+\alpha^{1}_{by}\cdot\sigma_{by}$
(12)
$\sigma_{2}=-\alpha^{2}_{br}\cdot\sigma_{br}+\alpha^{2}_{by}\cdot\sigma_{by}$
(13)

Here, for example, the $\alpha^{1}_{br}$ denotes for the stress concentration
factor at point 1 under bearing. Bearing stress $\sigma_{br}=P_{br}/dt$ and
by-pass stress $\sigma_{by}=P_{by}/(w-d)t$. As a result of the symmetric
position of point 1 and point 2, $\alpha^{1}_{by}$ equals to
$\alpha^{2}_{by}$. Therefore, bolt stress can be calculated from Eq.12 and
Eq.13 as,

$\sigma_{bolt1}=\sigma_{br}=\frac{\sigma_{1}-\sigma_{2}}{\alpha^{1}_{br}+\alpha^{2}_{br}}$
(14)

For all the three bolts, similar equations can be obtained, thus making the
bolt load ratio as,

$ratio1=\frac{\sigma_{bolt1}}{\sigma_{bolt1}+\sigma_{bolt2}+\sigma_{bolt3}}=\frac{\sigma_{1}-\sigma_{2}}{\sigma_{1}-\sigma_{4}}$
(15)
$ratio2=\frac{\sigma_{bolt2}}{\sigma_{bolt1}+\sigma_{bolt21}+\sigma_{bolt3}}=\frac{\sigma_{2}-\sigma_{3}}{\sigma_{1}-\sigma_{4}}$
(16)
$ratio2=\frac{\sigma_{bolt3}}{\sigma_{bolt1}+\sigma_{bolt21}+\sigma_{bolt3}}=\frac{\sigma_{3}-\sigma_{4}}{\sigma_{1}-\sigma_{4}}$
(17)

It can be seen from our method that the requirement for these four stain
gauges is that they are placed symmetric to each other at the corresponding
bolt hole. In our experiments, all the strain gauges are 30mm distant away
from hole center. The bolt hole clearances are adjusted by a careful
measurement of bolt diameter and the hole dimension cut by waterjet. The
tightening torques of three bolts are ensured by a ETC Series digital torque
wrench from ADEMA company. The fabricated three-bolt joints and the
experimental set-up are shown in Fig.14. In the experiments, a total external
load of 30kN was applied and kept when the strain values at the four strain
gauges at each joined plate were recorded and then calculated according to
Eq.15-Eq.17.

A total of four multi-bolt joints were made and follows two different
parameters as (0.2mm, 0.2mm, 0.2mm, 3.5Nm, 3.5Nm, 3.5Nm) and (1.0mm, 0.9mm,
1.0mm, 3.2Nm, 10.0Nm, 1.2Nm). The experimental results calculated using
Eq.15-Eq.17 and results by circuit model are shown in Fig.15. It can be shown
that for these two configurations, results by circuit model agree well with
the experimental method, especially at the middle bolt 2. For the bolt 1 and
bolt 3, relatively large error is shown. This is due to the secondary bending
effect as shown in Fig.14. Normally, this bending would have effect on value
of strain gauges 1 and 4. But at the strain gauge 2 and 3, this bending is
restrained by the contact brought by tightening torques. Therefore, a better
accuracy is shown in for the mid bolt. Comparing the results of these two
configurations, it can be found that for the first group of specimens, bolt 1
and bolt 3 carry clearly more load than bolt 2. With parameters (1.0mm, 0.9mm,
1.0mm, 3.2Nm, 10.0Nm, 1.2Nm) provided by optimization method ”GA+Circuit”,
results of the second group of specimens give a better even bolt-load-
distribution as bolt 2 carries more load. This is achieved by relatively
smaller the clearance and larger value of the tightening torque at bolt 2.
However, due to the complex multi-stage pattern brought by friction effect,
the mechanism of the combination of clearances and torques is not
straightforward and this is where our machine learning model can play an
important role.

Figure 13: Strain gauges placed on the multi-bolt joint structures and the
stress analysis Figure 14: (a)The three-bolt joints cut by waterjet machine
and assembled by a digital torque wrench (b)Experimental set-up (c)Secondary
bending effect Figure 15: Experimental validation for the optimized design at
an external tensile load of 30kN: comparison between results by circuit model
and experiments under a three-bolt joint with (a)uniform parameters (b)
parameters provided by optimization framework

##  5 Conclusion

To conclude, in this paper, we have proposed an effective optimization method
for the design of multi-bolt joint through a combination of a simple machine
learning model and a novel equivalent circuit model.

The idea of this equivalent circuit model for determining load distribution in
multi-bolt joint comes from the similarity between circuit system and linear
mechanic system. According to this similarity, the load shared by each part of
the joint structure is equal to the current flow in circuit model. Apart from
the advantage that circuit model shows better computational efficiency, effect
of frictions can also be reflected directly in this model by choosing various
parameters indicating bolt-hole clearances and tightening-torques. Moreover,
this circuit model is not limited to the single-bolt or three-bolt joint shown
in the paper. It can be easily extended to joint with more bolts through
serial and parallel modules.

While the comparative study between identical and random parameters of
clearances and torques has clearly shown that this friction effect has a great
influence on the load distribution ratio on each bolt, we are focusing on
reaching the minimum bolt load unevenness by carefully alter these two
significant parameters of each bolt.

The machine-learning based framework here aims to provide both the forward
response predictor and the backward problem of optimization. First, it is used
to generate a symbol of load distribution unevenness from given set of
parameters. Since that the dimensions of input and output is small, a very
simple machine learning model using neural network with low computational cost
is successfully established and validated. We then form a database which is
composed of three bolt-hole clearances in the range (0,2mm), three tightening
torques in the range (0,15Nm), and the corresponding symbol of unevenness of
bolt-load-distribution. Then, with these database, for the optimization
problem of load distribution unevenness, machine learning model could provide
optimal parameters for achieving minimum unevenness. In addition, by a
collaboration with optimization algorithms like GA and PSO, our results show
that the computational efficiency is greatly increased, which is the most
impressive advantages of our framework.

An experimental method of simply capturing the bolt load distribution using
four strain gauges is introduced. This is based on an assumption of linear
relationship between bolt load and stress at a certain position. By adopting
this experimental method, the output of our optimization framework is
validated.

The machine learning based framework can be a general procedure for
accelerating the inverse problem elucidating structure-property relationship.
A possible application could be problems where a complicated relationship is
hidden between input and output, or a vast design space is involved.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Acknowledgement

The authors are grateful to the support from NSFC/RGC Joint Research Scheme of
Hong Kong (Grant #: N_HKUST 631/18), Nanhai-HKUST Program (Grant #:
FSNH-18FYTRI01).

## Author contribution

Jinglei Yang and Zhidong Guan proposed and designed the project. Cheng Qiu
provided the idea of circuit model and performed benchmark calculation.
Fengyang Jiang, Yuzi Han, Cheng Qiu developed machine-learning model. Shanyi
du supervised the project. All authors contributed to the analysis of data as
well as writing and revising of the manuscript.

## References

  * [1] Srinivasa D Thoppul, Joana Finegan, and Ronald F Gibson.  Mechanics of mechanically fastened joints in polymer–matrix composite structures–a review.  Composites Science and Technology, 69(3-4):301–329, 2009.
  * [2] K B Katnam, LFM Da Silva, and TM Young.  Bonded repair of composite aircraft structures: A review of scientific challenges and opportunities.  Progress in Aerospace Sciences, 61:26–42, 2013.
  * [3] Karthik Prasanna Raju, Kobye Bodjona, Gyu-Hyeong Lim, and Larry Lessard.  Improving load sharing in hybrid bonded/bolted composite joints using an interference-fit bolt.  Composite Structures, 149:329–338, 2016.
  * [4] Pedro Ponces Camanho and Michel Lambert.  A design methodology for mechanically fastened joints in laminated composite materials.  Composites Science and Technology, 66(15):3004–3020, 2006.
  * [5] Yi Xiao and Takashi Ishikawa.  Bearing strength and failure behavior of bolted composite joints (part i: Experimental investigation).  Composites Science and Technology, 65(7-8):1022–1031, 2005.
  * [6] Yi Xiao and Takashi Ishikawa.  Bearing strength and failure behavior of bolted composite joints (part ii: modeling and simulation).  Composites science and technology, 65(7-8):1032–1043, 2005.
  * [7] Alfred Puck and Helmut Schürmann.  Failure analysis of frp laminates by means of physically based phenomenological models.  Composites science and technology, 62(12-13):1633–1662, 2002.
  * [8] ST Pinho, L Iannucci, and P Robinson.  Physically based failure models and criteria for laminated fibre-reinforced composites with emphasis on fibre kinking. part ii: Fe implementation.  Composites Part A: Applied Science and Manufacturing, 37(5):766–777, 2006.
  * [9] MA McCarthy, CT McCarthy, and GS Padhi.  A simple method for determining the effects of bolt–hole clearance on load distribution in single-column multi-bolt composite joints.  Composite Structures, 73(1):78–87, 2006.
  * [10] J Lecomte, C Bois, H Wargnier, and J-C Wahl.  An analytical model for the prediction of load distribution in multi-bolt composite joints including hole-location errors.  Composite Structures, 117:354–361, 2014.
  * [11] Yuxing Yang, Xueshu Liu, Yi-Qi Wang, Hang Gao, Yongjie Bao, and Rupeng Li.  An enhanced spring-mass model for stiffness prediction in single-lap composite joints with considering assembly gap and gap shimming.  Composite Structures, 187:18–26, 2018.
  * [12] Fengrui Liu, Xuheng Lu, Libin Zhao, Jianyu Zhang, Ning Hu, and Jifeng Xu.  An interpretation of the load distributions in highly torqued single-lap composite bolted joints with bolt-hole clearances.  Composites Part B: Engineering, 138:194–205, 2018.
  * [13] Libin Zhao, Ziang Fang, Fengrui Liu, Meijuan Shan, and Jianyu Zhang.  A modified stiffness method considering effects of hole tensile deformation on bolt load distribution in multi-bolt composite joints.  Composites Part B: Engineering, 171:264–271, 2019.
  * [14] CT McCarthy and PJ Gray.  An analytical model for the prediction of load distribution in highly torqued multi-bolt composite joints.  Composite Structures, 93(2):287–298, 2011.
  * [15] Cheng Qiu, Zhidong Guan, Shanyi Du, and Zengshan Li.  An improved model for predicting stiffness of single-lap composites bolted joints using matlab/simulink.  Mechanics of Advanced Materials and Structures, pages 1–11, 2019\.
  * [16] CT McCarthy and MA McCarthy.  Three-dimensional finite element analysis of single-bolt, single-lap composite bolted joints: Part ii—-effects of bolt-hole clearance.  Composite Structures, 71(2):159–175, 2005.
  * [17] Jonathan Schmidt, Mário RG Marques, Silvana Botti, and Miguel AL Marques.  Recent advances and applications of machine learning in solid-state materials science.  npj Computational Materials, 5(1):1–36, 2019.
  * [18] Chun-Teh Chen and Grace X Gu.  Machine learning for composite materials.  MRS Communications, 9(2):556–566, 2019.
  * [19] Rampi Ramprasad, Rohit Batra, Ghanshyam Pilania, Arun Mannodi-Kanakkithodi, and Chiho Kim.  Machine learning in materials informatics: recent applications and prospects.  npj Computational Materials, 3(1):1–13, 2017.
  * [20] Navid Zobeiry, Johannes Reiner, and Reza Vaziri.  Theory-guided machine learning for damage characterization of composites.  Composite Structures, page 112407, 2020.
  * [21] Logan Ward, Ankit Agrawal, Alok Choudhary, and Christopher Wolverton.  A general-purpose machine learning framework for predicting properties of inorganic materials.  npj Computational Materials, 2(1):1–7, 2016.
  * [22] Liam Huber, Raheleh Hadian, Blazej Grabowski, and Jörg Neugebauer.  A machine learning approach to model solute grain boundary segregation.  npj Computational Materials, 4(1):1–8, 2018.
  * [23] Zijiang Yang, Yuksel C Yabansu, Reda Al-Bahrani, Wei-keng Liao, Alok N Choudhary, Surya R Kalidindi, and Ankit Agrawal.  Deep learning approaches for mining structure-property linkages in high contrast composites from simulation datasets.  Computational Materials Science, 151:278–287, 2018.
  * [24] Chunping Ma, Zhiwei Zhang, Benjamin Luce, Simon Pusateri, Binglin Xie, Mohammad H Rafiei, and Nan Hu.  Accelerated design and characterization of non-uniform cellular materials via a machine-learning based framework.  npj Computational Materials, 6(1):1–8, 2020.
  * [25] Ruoqian Liu, Abhishek Kumar, Zhengzhang Chen, Ankit Agrawal, Veera Sundararaghavan, and Alok Choudhary.  A predictive machine learning approach for microstructure optimization and materials design.  Scientific reports, 5(1):1–12, 2015.
  * [26] Chun-Teh Chen and Grace X Gu.  Generative deep neural networks for inverse materials design using backpropagation and active learning.  Advanced Science, 7(5):1902607, 2020.
  * [27] Grace X Gu, Chun-Teh Chen, Deon J Richmond, and Markus J Buehler.  Bioinspired hierarchical composite design using machine learning: simulation, additive manufacturing, and experiment.  Materials Horizons, 5(5):939–945, 2018.
  * [28] P Ramasamy and S Sampathkumar.  Prediction of impact damage tolerance of drop impacted wgfrp composite by artificial neural network using acoustic emission parameters.  Composites Part B: Engineering, 60:457–462, 2014.
  * [29] Anh Tran, John A Mitchell, Laura Swiler, and Tim Wildey.  An active learning high-throughput microstructure calibration framework for solving inverse structure-process problems in materials informatics.  Acta Materialia, 2020.
  * [30] CT McCarthy, MA McCarthy, WF Stanley, and VP Lawlor.  Experiences with modeling friction in composite bolted joints.  Journal of composite materials, 39(21):1881–1908, 2005.
  * [31] Inc Abaqus.  Abaqus user’s manual.  ABAQUS Ine, Dassault Systemes, Rhode Island, USA, 2007.
  * [32] Fengrui Liu, Libin Zhao, Saqib Mehmood, Jianyu Zhang, and Binjun Fei.  A modified failure envelope method for failure prediction of multi-bolt composite joints.  Composites Science and Technology, 83:54–63, 2013.


# Galactic Clustering Under Power-law Modified Newtonian Potential

Abdul W. Khanday<EMAIL_ADDRESS>Department of Physics, National
Institute of Technology, Srinagar, Kashmir-190006, India  Sudhaker Upadhyay
<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of
Physics, K.L.S. College, Nawada, Bihar 805110, India  Department of Physics,
Magadh University, Bodh Gaya, Bihar 824234, India  Inter-University Centre for
Astronomy and Astrophysics (IUCAA) Pune, Maharashtra 411007, India  School of
Physics, Damghan University, P. O. Box 3671641167, Damghan, Iran  Prince A.
Ganai<EMAIL_ADDRESS>Department of Physics, National Institute of
Technology, Srinagar, Kashmir-190006, India

###### Abstract

We estimate galaxy clustering under a modified gravitational potential. In
particular, the modifications in gravitational potential energy occur due to a
power-law and cosmological constant terms. We derive a canonical partition
function for the system of galaxies interacting under such a modified
gravitational potential. Moreover, we compute various thermodynamical equation
of states for the system. We do comparative analysis in order to emphasize the
effect of corrections on thermodynamics of the system. Interestingly, the
modifications in thermodynamical quantities are embedded in clustering
parameter only.

##  1 Introduction

Since last decade, the substantial progresses have been made in the
understandings of galaxy clusters, from their internal structure and evolution
to their place in the large scale structure of the universe. All these
progresses are due to the stupendous improvements in the theoretical modeling
and numerical simulation, viz-a-viz abundance of new information provided by
multi-wavelength surveys of the universe. Various theories of cosmological
many-body distribution function have been developed from the thermodynamic
point of view. Benkestien [1], Hawking [2] and Unruh [3] originated the
relation between relativity and thermodynamics. Later, Jacobson [4] introduced
the Einstein equation as thermodynamic equation of state. In fact, the
Einstein equation is derived from the proportionality of entropy and the
horizon area together with the fundamental relation $\delta Q=TdS$, where
$\delta Q$ and $T$ are interpreted as the energy flux and Unruh temperature.
Recently Verlinde [5] proposed an entropic origin of gravity and interpreted
gravity as an entropic force. In this regard, it is argued that the central
notion needed to derive gravity is information. Employing the holographic
principle and equipartition law of energy, Newton’s law of gravitation,
Poisson’s equation and Einstein’s field equations are successfully obtained.
Verlinde had used area-entropy relation of black holes in Einstein’s gravity
(i.e. $S=\frac{A}{4l_{p}^{2}}$) to get the Newtons law, where $S,A$ and
$l_{P}$ refer the entropy of the black hole, the area of horizon and the
Planck length, respectively.

The study of modified gravity on the cosmological scales is an active area of
present research. Recently, the corrected gravitational potential plays a
vital role in estimation of the total mass of a sample of 12 clusters of
galaxies which provides a better fit to the mass of visible matter [6]. At
large distances, the modification in Newtonian potential occurs due to the
propagation of gravity into the bulk [7]. Recently, Sheykhi and Hendi studied
the effect of power-law corrections in entropy to the Newton’s law [8]. These
corrections appear due to entanglement of quantum field inside and outside of
horizon [9]. It is also found that a viable source of black hole entropy is
quantum entanglement of degrees of freedom inside and outside the horizon.
Also, the black hole entropy is found directly proportional to the surface
area of the sphere when the field is in the ground state. But when the field
is in a superposition of ground and excited states, a correction term
proportional to a fractional power of area appears. These corrections are
negligible for large horizon areas. Another quantum correction to Newton’s
law, a logarithmic correction, has also been studied which appears due to the
result of thermal equilibrium fluctuations and the quantum fluctuations of
loop quantum gravity [10, 11, 12]. Recently, the effect of this logarithmic
correction on the galaxy clustering has also been studied [13]. The clustering
of galaxies has been studied under the modified gravity under various
potentials [14, 15, 16, 17, 18, 19, 20, 21].

The power-law corrected entropy leads to modification in Newton’s
gravitational potential by adopting the viewpoint of gravity as an entropic
force. We consider a power-law corrected entropy of the following form [9,
22]:

$S=\dfrac{A}{4l_{p}^{2}}[1-K_{\alpha}A^{1-\frac{\alpha}{2}}],$ (1)

where $\alpha$ is a dimensionless constant whose value is not confirmed yet
and the parameter $K_{\alpha}$ that depends on the power $\alpha$ of the
entropy correction as following:

$K_{\alpha}=\frac{\alpha(4\pi)^{\frac{\alpha}{2}-1}}{(4-\alpha)r_{c}^{2-\alpha}},$

here $r_{c}$ denotes crossover scale. The Boltzmann constant is set unit here.
The last term of equation (1) appears due to the superposition of ground state
and excited state wave function of the field. Interestingly, Verlinde [5]
proposed an entropic origin of gravity and interpreted gravity as an entropic
force. The gravity derived here was from the notion of information associated
with matter and its location measured in terms of entropy. The Newton’s law of
gravitational force $F$ is related to the entropy $S$ of the system as [8]

$\displaystyle F=-4l_{p}^{2}\frac{GM^{2}}{r^{2}}\frac{\partial S}{\partial
A}.$ (2)

For the power-law corrected entropy given in (1), this relation leads to the
following modified Newton’s force [8]:

$F=-\frac{GM^{2}}{r^{2}}\left[1-\frac{\alpha}{2}\left(\frac{r_{c}}{r}\right)^{\alpha-2}\right],$
(3)

which coincides with the original Newtons law when $\alpha$ is set to zero.
Keeping in view the attractive nature of gravity, we should have $F<0$. This
requires,

$\alpha<2\left(\frac{r}{r_{c}}\right)^{\alpha-2}.$

In this paper, we study the effect of power-law corrected Newtonian potential
on the clustering of galaxies. We also consider the effect of dark matter on
the formation of galaxy clusters via the incorporation of cosmological
constant in the Newtonian potential. This is because of the role played by
cosmological constant $\Lambda$ in the expansion of universe [5]. In order to
study the effects of all the corrections made in gravitational potential in
galaxy clustering, we first derive the $N$-body partition function by
evaluating configuration integrals recursively. The resulting partition
function is employed to calculate the various thermodynamic quantities viz the
Helmholtz free energy, entropy, pressure, internal energy, and chemical
potential which possess deviations from their original values due to the
incorporation of corrections. Remarkably, a modified clustering parameter
emerges naturally from the corrected equations of state. Furthermore, we
derive the probability distribution function assuming that the system is in a
quasi equilibrium state as described by the grand canonical ensemble. The
expression of distribution function embeds modified clustering parameter. A
comparative analysis is made with the original distribution function to study
the deviation due to corrections.

The paper is organized as following. In section 2, we consider a power-law and
cosmological constant modified gravitational potential to calculated partition
function of the system of galaxies and galaxy clusters. With the help of
resulting partition function, the various thermodynamical equation of states
are calculated in section 3. The study the distribution of galaxies and
galaxy’s clusters under the modified gravitational potential, we estimate
distribution function in section 4. Finally, we draw concluding remarks in the
last section.

##  2 Interaction of galaxy clusters under modified potential

In this section, we derive a modified Newtonian potential due to power-law
corrected force and estimate the partition function.

###  2.1 A modified gravitation potential

Utilized standard definition $\Phi=-\int Fdr$, we calculate a power-law
modified gravitational potential corresponding to power-law corrected Newton’s
law (3) as

$\Phi=-GM^{2}\left[\frac{1}{r}+\frac{\alpha}{2\left(\alpha+1\right)}r_{c}^{\alpha-2}r^{\alpha+1}\right].$
(4)

This potential can arise in modified gravity theories like $f(R)$ gravity. In
fact, the power-law entropy corrected Friedmann equation is derived using the
first-law on the apparent horizon [8].

To get more realistic results, one can not ignore the cosmological constant
term $-\frac{1}{6}\Lambda r$ [23] in potential at the cosmological scale as
cosmological constant is responsible for the expansion of the Universe through
a repulsive force. Therefore, the exact gravitation potential is given by

$\Phi=-GM^{2}\left[\frac{1}{r}+\frac{\alpha}{2\left(\alpha+1\right)}r_{c}^{\alpha-2}r^{\alpha+1}+\frac{1}{6}\frac{\Lambda
r^{2}}{GM^{2}}\right].$ (5)

This is a final gravitational potential where second and third terms
correspond to the power-law and cosmological constant correction terms
respectively.

###  2.2 Generating functional of galaxies cluster under modified gravity

Next, we estimate the partition function for the system of galaxies under the
modified gravitational potential. Here it is assumed that system of galaxies
follows a statistically homogeneous distribution over large regions, which
consists of an ensemble of cells with equal volume $V$ and equal average
density. Let us begin by writing the general partition function for the system
comprised with $N$ galaxies of equal mass $M$, momenta $P_{i}$ and average
temperature $T$ as

$Z_{N}\left(T,V\right)=\dfrac{1}{\lambda^{3N}N!}\int d^{3N}Pd^{3N}r\
\exp\left[-\frac{1}{T}\left(\sum_{i=1}^{N}\frac{P_{i}^{2}}{2M}+\Phi\left(r_{1},r_{2},...r_{N}\right)\right)\right],$
(6)

where $N!$ appears due to distinguishable nature of galaxies and $\lambda$ is
a normalization constant for the phase space volume cell. Upon integration
over momentum space, this further simplifies to

$Z_{N}\left(T,V\right)=\frac{1}{N!}\left(\frac{2\pi
MT}{\lambda^{2}}\right)^{3N/2}Q_{N}\left(T,V\right),$ (7)

where the configuration integral, $Q_{N}\left(T,V\right)$, has the following
form:

$Q_{N}\left(T,V\right)=\int...\int\prod_{1\leq i<j\leq
N}\left(1+\frac{\Phi}{T}\right)d^{3N}r,$ (8)

Here we neglected the higher-order terms of potential as the system of
galaxies is still clustering.

Here we note that for the point-mass galaxies (i.e., for $r=0$), the potential
energy and consequently partition function diverges. In order to remove this
discrepancy, we assume galaxies of extended nature (i.e., galaxies with
halos). For this we introduce a softening parameter $\epsilon$, which assures
the finite size of galaxies. The value of this softening parameter ranges
$0.01\leq\epsilon\leq 0.05$. In order to estimate partition function, the
modified potential incorporates softening parameter appropriately as follows,

$\Phi\left(\epsilon\right)=-GM^{2}\left[\frac{1}{\left(r^{2}+\epsilon^{2}\right)^{\frac{1}{2}}}+\frac{\alpha
r_{c}^{\alpha-2}}{2\left(1-\alpha\right)}r^{1-\alpha}+\frac{1}{6}\frac{\Lambda
r^{2}}{GM^{2}}\right].$ (9)

The second and third terms do not require softening parameter to be introduces
as corresponding potentials do not diverge.

Now, we estimate configuration integral $Q_{N}\left(T,V\right)$ iteratively.
The configuration integral for a single (spherically) galaxy of radius $R$ is
given by

$Q_{1}\left(T,V\right)=V.$ (10)

For system of two galaxies, the configuration integral is given by

$Q_{2}\left(T,V\right)=4\pi V\int dR\
R^{2}\left[1+\frac{GM^{2}}{T}\left(\frac{1}{\left(R^{2}+\epsilon^{2}\right)^{\frac{1}{2}}}+\frac{\alpha
r_{c}^{\alpha-2}}{2\left(1-\alpha\right)}R^{1-\alpha}+\frac{1}{6}\frac{\Lambda
R^{2}}{GM^{2}}\right)\right].$ (11)

Here, the double integral reduces to a single integral by fixing the position
of one galaxy. This result can also be obtained by considering the fact that
expansion of universe exactly cancels the effect of the long-range mean
gravitational field on the particle motions [24]. Eq. (11) further simplifies
to

$\displaystyle Q_{2}\left(T,V\right)$ $\displaystyle=$ $\displaystyle
V^{2}\left[1+\frac{3}{2}\frac{GM^{2}}{RT}\left(\sqrt{1+\frac{\epsilon^{2}}{R^{2}}}+\frac{\epsilon^{2}}{R^{2}}\ln\frac{{\epsilon}/{R}}{1+\sqrt{1+\frac{\epsilon^{2}}{R^{2}}}}\right.\right.$
(12) $\displaystyle+$
$\displaystyle\left.\left.\frac{3\alpha}{2\left(1-\alpha\right)\left(4-\alpha\right)}r_{c}^{\alpha-2}R^{1-\alpha}+\frac{\Lambda
R^{3}}{15GM^{2}}\right)\right].$

More compactly, this can be written as

$Q_{2}\left(T,V\right)=V^{2}\left[1+\frac{3}{2}\left(\zeta+\gamma+\beta\right)\frac{GM^{2}}{RT}\right],$
(13)

where

$\displaystyle\zeta$
$\displaystyle=\sqrt{1+\frac{\epsilon^{2}}{R^{2}}}+\frac{\epsilon^{2}}{R^{2}}\ln\frac{{\epsilon}/{R}}{1+\sqrt{1+\frac{\epsilon^{2}}{R^{2}}}},$
$\displaystyle\gamma$
$\displaystyle=\frac{3\alpha}{2\left(1-\alpha\right)\left(4-\alpha\right)}r_{c}^{\alpha-2}R^{1-\alpha},$
$\displaystyle\beta$ $\displaystyle=\frac{\Lambda R^{3}}{15GM^{2}}.$

Here, we note that $\alpha$ in the $\gamma$ takes the value $0.624<\alpha<2$
because within this limit only a phantom accelerating universe can be derived
which is compatible with the observations [25].

Next, we scale the temperature $T$ and radius $R$ as $T\rightarrow\eta^{-1}T$
and $R\rightarrow\eta R$, which leads to the dimensionless factor
$\frac{GM^{2}}{RT}$ scale invariant. Therefore, we can scale
$\frac{GM^{2}}{RT}\rightarrow\left(\frac{GM^{2}}{RT}\right)^{3}=\frac{3}{2}\left(\frac{GM^{2}}{T}\right)^{3}\bar{\rho}:=x$
[26]. Thus, Eq. (13) finally reduces to

$Q_{2}\left(T,V\right)=V^{2}\left[1+\left(\zeta+\gamma+\beta\right)x\right].$
(14)

Following this procedure iteratively, the configuration integral for $N$
galaxies is obtained as

$Q_{N}\left(T,V\right)=V^{N}\left[1+\left(\zeta+\gamma+\beta\right)x\right]^{N-1}.$
(15)

By inserting this value of configuration integral of $N$ galaxies into the the
partition function (7), we achieve the expression of partition function for
$N$ galaxies gravitating under modified gravity as

$Z_{N}\left(T,V\right)=\frac{1}{N!}\left(\frac{2\pi
MT}{\lambda^{2}}\right)V^{N}\left[1+\left(\zeta+\gamma+\beta\right)x\right]^{N-1}.$
(16)

##  3 Thermodynamical equations of state

In this section, we derive various more exact equations of states for the
system of galaxies interacting under modified potential. More precisely, we
derive the Helmholtz free energy, entropy, internal energy, pressure and
chemical potential. We also emphasize the effect of corrections in Newton’s
potential on these equation of states.

###  3.1 Helmholtz free energy

The Helmholtz free energy can be calculated from the partition function using
following definition:

$F=-T\ln Z_{N}(T,V).$ (17)

Therefore, the Helmholtz free energy for our system of galaxies is calculated
by

$\displaystyle F=NT\left(\ln\frac{N}{V}T^{-3/2}\right)-NT-
NT\ln\left[1+\left(\zeta+\gamma+\beta\right)x\right]-\frac{3}{2}NT\ln\left(\frac{2\pi
M}{\lambda^{2}}\right).$ (18)

Here we have made following assumption $(N-1)\approx N$, since $N$ is very
large. This equation can also be written as

$F=NT\left(\ln\frac{N}{V}T^{-3/2}\right)-NT+NT\ln\left[1-\frac{(\zeta+\gamma+\beta)x}{1+(\zeta+\gamma+\beta)x}\right]-\frac{3}{2}NT\ln\left(\frac{2\pi
M}{\lambda^{2}}\right).$ (19)

Figure 1: Helmholtz free energy vs. number of galaxies. Left: red, blue, and
green lines correspond to $\gamma x=0,0.5$ and $1$, respectively, with $\zeta
x=1$, $\beta x=0$. Right: red, blue, and green lines correspond to $\gamma
x=0,0.5$ and $1$, respectively, with $\zeta x=1$, $\beta x=1$. Rest of the
parameters are unit here.

Figure 2: Behavior of Helmholtz free energy versus particle number for
different values of $M$. We set all the parameters along with $\zeta x$,
$\beta x$ and $\gamma x$ to unit. Here, red, blue, and green lines correspond
to $M=1$, $M=5$ and $M=10$, respectively.

In order to study the behavior of Helmholtz free energy with respect to
particle number, we plot a graph 1 and 2. From the figure, we see that
Helmholtz free energy is a decreasing function of $N$. It means as long as
number of galaxies increases the Helmholtz free energy decreases. It is also
evident that the corrections due to power-law and cosmological constant terms
make the free energy more negative. Also, the free energy decreases with
increase in the mass of galaxies.

By comparing this equation of state to its original form given in Ref. [24],
we conclude that the corrections in Newton’s potential are apparent in
clustering parameter only. The new clustering parameter corresponding to
power-law and cosmological constant corrections is obtained as

$b_{\star}=\frac{(\zeta+\gamma+\beta)x}{1+(\zeta+\gamma+\beta)x}.$ (20)

This modified clustering parameter can be expressed in terms of original
clustering parameter $b_{\epsilon}=\frac{\zeta x}{1+\zeta x}$ [24] as follows

$b_{\star}=\frac{b_{\epsilon}(1-\gamma x-\beta
x)+(\gamma+\beta)x}{1+(\gamma+\beta)x-b_{\epsilon}(\gamma+\beta)x}.$ (21)

###  3.2 Entropy

Let us calculate entropy of the system which a very important thermodynamical
quantity. For a given Helmholtz free energy (19), the entropy can easily be
calculated as

$\displaystyle
S=-N\ln\left(\frac{N}{V}T^{-3/2}\right)-N\ln[1-b_{\star}]-3Nb_{\star}+\frac{5}{2}N+\frac{3}{2}N\ln\frac{2\pi
M}{\lambda^{2}}.$ (22)

Here we have utilized standard entropy definition $S=-\left(\frac{\partial
F}{\partial T}\right)_{N,V}$. The specific entropy (entropy per galaxy) is
given by

$\frac{S}{N}=-\ln\left(\frac{N}{V}T^{-3/2}\right)-\ln[1-b_{\star}]-3b_{\star}+S_{0},$
(23)

where $S_{0}=\frac{5}{2}+\frac{3}{2}\ln\frac{2\pi M}{\lambda^{2}}$ is a
constant.

Figure 3: Specific entropy ($S/N$) versus particle number ($N$) with (right)
and without (left) dark energy contributions. Left: red, blue, and green lines
correspond to $\gamma x=0,0.5$ and $1$, respectively, with $\zeta x=1$ and
$\beta x=0$. Right: red, blue, and green lines correspond to $\gamma x=0,0.5$
and $1$, respectively with $\zeta x=1$ and $\beta x=1$. Rest of the parameters
are unit here.

Figure 4: Behavior of specific entropy ($S/N$) versus particle number ($N$)
for different values of $M$. We set all the parameters along with $\zeta x$,
$\beta x$ and $\gamma x$ to unit. Here, red, blue, and green lines correspond
to $M=1$, $M=5$ and $M=10$, respectively.

The behavior of specific entropy with respect to number of galaxies can be
seen in Figs. 3 and 4. The specific entropy is also a positive valued
decreasing function with number of galaxies. The power-law corrected term
makes specific entropy smaller in absence of dark energy. However, the effect
of power-law correction is not significant in the presence of dark energy but
makes specific entropy bit positive. The specific entropy increases
significantly along with increase in mass of galaxies.

###  3.3 Internal Energy

The internal energy of the system is defined as $U=F+TS$. For a given
Helmholtz free energy (19) and entropy (22), the internal energy is calculated
by

$U=\frac{3}{2}NT\left[1-2\frac{\left(\zeta+\gamma+\beta\right)x}{1+\left(\zeta+\gamma+\beta\right)x}\right]=\frac{3}{2}NT\left[1-2b_{\star}\right].$
(24)

Figure 5: Internal energy ($U$) versus particle number ($N$) with (right) and
without (left) dark energy contributions. Left: red, blue, and green lines
correspond to $\gamma x=0,0.5$ and $1$, respectively, with $\zeta x=2$ and
$\beta x=0$. Right: red, blue, and green lines correspond to $\gamma x=0,0.5$
and $1$, respectively with $\zeta x=2$ and $\beta x=1$. Rest of the parameters
are unit here.

From the figure 5, we observe that the internal energy of the system which is
negative valued becomes more negative under the effect of correction in
gravitational potential.

###  3.4 Pressure

The standard definition of pressure is given by $P=-\left(\frac{\partial
F}{\partial V}\right)_{N,T}.$ Exploiting expression (19), the pressure for a
gravitating system under modified gravity is given by

$P=\frac{NT}{V}\left[1-b_{\star}\right].$ (25)

Figure 6: Pressure ($P$) versus particle number ($N$) with (right) and without
(left) dark energy contributions. Left: red, blue, and green lines correspond
to $\gamma x=0,0.5$ and $1$, respectively, with $\zeta x=2$ and $\beta x=0$.
Right: red, blue, and green lines correspond to $\gamma x=0,0.5$ and $1$,
respectively with $\zeta x=2$ and $\beta x=1$. Rest of the parameters are unit
here.

The pressure is a linear function of number of the galaxies, which means that
pressure increases as long as number of galaxies increases. It is obvious from
the figure 6 that the presence of correction terms negate the pressure of the
system.

###  3.5 Chemical Potential

The chemical potential which measure exchange of galaxies can be calculated
from the formula $\mu=-\left(\frac{\partial F}{\partial N}\right)_{V,T}$. So,
it is a matter of calculation to derive chemical potential for this system as

$\mu=T\left(\ln\frac{N}{V}T^{-3/2}\right)-T\ln\left[1+(\zeta+\gamma+\beta)x\right]-T\frac{(\zeta+\gamma+\beta)x}{1+(\zeta+\gamma+\beta)x}-\frac{3}{2}T\ln\left(\frac{2\pi
M}{\lambda^{2}}\right).$ (26)

This further simplifies to

$\mu=T\left(\ln\frac{N}{V}T^{-3/2}\right)+T\ln\left[1-b_{\star}\right]-Tb_{\star}-\frac{3}{2}T\ln\left(\frac{2\pi
M}{\lambda^{2}}\right).$ (27)

Figure 7: Chemical potential ($\mu(N)$) versus particle number ($N$) with
(right) and without (left) dark energy contributions. Left: red, blue, and
green lines correspond to $\gamma x=0,0.5$ and $1$, respectively, with $\zeta
x=1$ and $\beta x=0$. Right: red, blue, and green lines correspond to $\gamma
x=0,0.5$ and $1$, respectively with $\zeta x=1$ and $\beta x=1$. Rest of the
parameters are unit here.

Figure 8: Behavior of chemical potential ($\mu(N)$) versus particle number
($N$) for different values of $M$. We set all the parameters along with $\zeta
x$, $\beta x$ and $\gamma x$ to unit. Here, red, blue, and green lines
correspond to $M=1$, $M=5$ and $M=10$, respectively.

The behavior of chemical potential with respect to number of galaxies can be
seen from figures 7 and 8. Interestingly, we observe that the chemical
potential is negative valued for the system of galaxies with small number of
galaxies. However, it becomes positive valued when number of galaxies
increases to a specific value. The presence of correction terms in clustering
parameter decreases the chemical potential of the system. The chemical
potential decreases with increase of the mass of galaxies.

At the end of this section, we conclude that the modified Newton’s potential
amounts changes to the clustering parameter $b_{\star}$ only, however the
basic structures of the equations remain intact.

##  4 General distribution function

In order to find the probability distribution function $F(N)$, which contains
void distribution as well as statistics of counts of the number of galaxies in
cells throughout the system. For the system of galaxies, wherein galaxies as
well as energy can cross the cell boundary, one has to estimate the grand
canonical partition function. The grand canonical partition function, a
weighted sum of all canonical partition functions, is defined by

$Z_{G}(T,V,z)=\sum_{N=0}^{\infty}e^{\frac{N\mu}{T}}Z_{N}(T,V),$ (28)

where $z$ is the fugacity. The grand partition function for the system of
galaxies is expressed in terms of thermodynamic variables as

$\ln Z_{G}=\frac{PV}{T}=N(1-b_{\star}).$ (29)

The probability distribution function $F(N)$ for finding $N$ galaxies in a
cell of volume $V$ and energy $U(N,V)$ is defined by

$F(N)=\frac{\sum_{i}e^{\frac{N\mu}{T}}e^{\frac{U_{i}}{T}}}{Z_{G}(T,V,z)}=\frac{e^{\frac{N\mu}{T}}Z_{N}(V,T)}{Z_{G}(T,V,z)}.$
(30)

Making use of Eqs. (16), (27) and (29), the distribution function for the
extended mass galaxies under modified potential is estimated as

$F(N)=\frac{\bar{N}}{N!}(1-b_{\star})[\bar{N}(1-b_{\star})+Nb_{\star}]^{N-1}e^{-Nb_{\star}-\bar{N}(1-b_{\star})}.$
(31)

This distribution function is structurally similar to those derived originally
by Saslaw and Hamilton [27] from thermodynamic point of view and by Ahmad and
Saslaw [24] from statistical point of view. Also the distribution function
derived from logarithmic and volume corrected Newtonian potential [13] has the
same general structure.

Figure 9: The distribution function $F(N)$ versus $N$ for $\bar{N}=10$ with
(right) and without (left) cosmological constant contributions. Left: red,
blue, and green lines correspond to $\gamma x=0,0.5$ and $1$, respectively,
with $\zeta x=1$ and $\beta x=0$. Right: red, blue, and green lines correspond
to $\gamma x=0,0.5$ and $1$, respectively with $\zeta x=1$ and $\beta x=1$.

The behavior of distribution function $F(N)$ versus $N$ can be seen from the
comparative analysis as given in Fig. 9. The presence of correction term
decreases the peak value of distribution function which occurs for system of
small number of galaxies. However, for the system of large number of galaxies
the correction terms increase the distribution function.

##  5 Discussion and conclusions

We have presented a study of galaxy clustering under the modified Newton’s
law. These modifications to Newton’s law incorporate a power-law entropic
corrections along with the inclusion of cosmological constant $\Lambda$ term,
in order to take into account the effect of quantum entanglement (a possible
source of black hole entropy) and dark energy respectively. Utilizing the
modified Newtonian potential, we have derived the corresponding canonical
partition function for the system of extended mass galaxies with the
assumption that the system is made of $N$ equal volume cells and average
particle density $\bar{\rho}$, which are statistically homogeneous over large
regions. In this regard, we have used of softening parameter $\epsilon$ in the
first term of Hamiltonian to get rid of the divergence of the Hamiltonian when
galaxies are considered as point like. This justifies that the extended nature
(finite size) of galaxies (or galaxies with halos). From the resulting
partition function, we have calculated various important thermodynamic
equations of state. Namely, these are free energy, entropy, internal energy,
pressure and chemical potential. The exact expressions of equations of state
contain a corrected correlation (clustering) parameter which emerges naturally
for the clusters of the galaxies with halos. The new clustering parameter
$b_{\star}$ reduces to the original parameter $b_{\epsilon}$ when $\gamma=0$
and $\beta=0$. Moreover, the distribution function is modified due to
correction in potential but has the similar structure as the original one.

## References

  * [1] J. D. Benkenstein, _Phys. Rev. D_ 7, 2333 (1973)
  * [2] S. W. Hawking, _Commun. Math. Phys._ 43, 199 (1975)
  * [3] W. G. Unruh, _Phys. Rev. D_ 14, 870, (1976)
  * [4] T. Jacobson, _Phys. Rev. Lett._ 75, 1260 (1995)
  * [5] E. P. Verlinde, _J. High Energy Phys._ 04, 029, (2011)
  * [6] S. Capozziello, E. De Filippis and V. Salzano, _MNRAS_ 394, 947 (2009)
  * [7] E. G. Floratos and G. K. Leontaris, _Phys. Lett. B_ 465, 95 (1999)
  * [8] A. Sheykhi and S. H. Hendi, _Phys. Rev. D_ 84, 044023 (2011)
  * [9] S. Das, S. Shankaranaraynan and S. Sur, _Phys. Rev. D_ 77, 064013 (2008)
  * [10] K. A. Meissner, _Classical Quantum Gravity_ 21, 5245 (2004)
  * [11] A. Ghosh and P. Mitra, _Phys. Rev. D_ 71, 027502 (2005)
  * [12] A. Chatterjee and Majumdar, _Phys. Rev. Lett._ 92, 141301 (2004)
  * [13] S. Upadhyay, _Phys. Rev. D_ 95, 043008 (2017)
  * [14] R. H. Sanders, _MNRAS_ 342, 901 (2003)
  * [15] M. Milgrom and R. H. Sanders, _Astroph. J._ 599, L25 (2003)
  * [16] M. Hameeda, S. Upadhyay, M. Faizal and A. F. Ali, _MNRAS_ 463, 3699 (2016)
  * [17] B. Pourhassan, S. Upadhyay, M. Hameeda and M. Faizal, _MNRAS_ 468, 3166 (2017)
  * [18] S. Upadhyay, S. Capozziello and B. Pourhassan, _Int. J. Mod. Phys. D_ 28, 1950027 (2019)
  * [19] S. Capozziello, M. Faizal, M. Hameeda, B. Pourhassan, V. Salzano and S. Upadhyay, _MNRAS_ 474, 2430 (2018)
  * [20] B. Pourhassan and S. Upadhyay, _Phys. Dark Univ._ 29, 100596 (2020)
  * [21] M. Hameeda, S. Upadhyay, M. Faizal, A. F. Ali and B. Pourhassan, _Phys. Dark Univ._ 19, 137 (2018)
  * [22] N. Radicella and D. Pavon, _Phys. Lett. B_ 691, 121 (2010)
  * [23] W. Rindler and M. Ishak, _Phys. Rev. D_ 76, 043006 (2007)
  * [24] F. Ahmad, W. C. Saslaw, and N. I. Bhat, _Astrophysics J._ 571, 576 (2002)
  * [25] K. Karami, A. Abdolmaleki, N. Sahraei and S. Ghaffari, _JHEP_ 1108, 150 (2011)
  * [26] L. D. Landau and E. M. Lifshtiz, _Statistical Physics_ (Pergamon Oxford, 1980)
  * [27] W. C. Saslaw and A. J. Hamilton, _Astrophysics J._ 13, 276 (1984)