# Modeling crack propagation in heterogeneous materials: Griffith’s law,
intrinsic crack resistance and avalanches

Subhadeep Roy<EMAIL_ADDRESS>Takahiro Hatano
<EMAIL_ADDRESS>Purusattam Ray<EMAIL_ADDRESS>PoreLab,
Department of Physics, Norwegian University of Science and Technology, N-7491
Trondheim, Norway.  Department of Earth and Space Science, Osaka University,
560-0043 Osaka, Japan.  The Institute of Mathematical Sciences, Taramani,
Chennai-600113, India.

###### Abstract

Various kinds of heterogeneity in solids including atomistic discreteness
affect the fracture strength as well as the failure dynamics remarkably. Here
we study the effects of an initial crack in a discrete model for fracture in
heterogeneous materials, known as the fiber bundle model. We find three
distinct regimes for fracture dynamics depending on the initial crack size. If
the initial crack is smaller than a certain value, it does not affect the
rupture dynamics and the critical stress. While for a larger initial crack,
the growth of the crack leads to breakdown of the entire system, and the
critical stress depends on the crack size in a power-law manner with a
nontrivial exponent. The exponent, as well as the limiting crack size, depend
on the strength of heterogeneity and the range of stress relaxation in the
system.

###### keywords:

Nucleation, Griffith’s criterion, Disordered systems, Fiber bundle model,
Avalanche statistics, Critical stress

##  1 Introduction

Determination of the strength of materials is a prime objective in material
science and applied mechanics. Traditionally the problem has theoretically
been addressed within the framework of continuum mechanics which relates the
motion of a crack to the applied loading [1]. For brittle solids, the problem
boils down to destabilization and subsequent growth of a pre-existing crack
under tension. Griffith (in 1920) [2] considered a single sharp crack in a
homogeneous elastic medium and suggested that the weakening of material by a
crack could be treated as an equilibrium problem in which the reduction of
strain energy, when the crack propagates, could be equated to the increase in
surface energy due to the increase in surface area. He found that the critical
stress $\sigma_{c}$ to cause a crack of length $l$, to extend is
$\sigma_{c}=(2Yg/\pi l)^{1/2}$, where $Y$ is Young’s modulus and $g$ is the
surface energy per unit area of the crack surface [2]. The Griffith criterion
has been extensively verified in engineering specimens containing cracks of
controlled length and is still used to estimate the surface energy in cracks
in brittle materials [1, 3].

However, the Griffith criterion is often found to lead to significant errors
in the estimation of surface energies. Even for perfectly brittle materials,
modifications of Griffith’s theory are needed to take into account the
discrete atomistic nature of the interactions [3]. The discreteness and
heterogeneities of materials lead to energy barriers at the crack tip which
can arrest crack motion: a phenomenon known as lattice trapping or intrinsic
crack resistance [4, 5, 6, 7, 8, 9]. Lattice trapping in crystalline materials
and intrinsic crack resistance in heterogeneous materials have received lots
of attention in the context of nucleation and growth of cracks in these
materials. One consequence of lattice trapping is the appearance of two length
scales in the systems [6, 4, 5, 10]: a large length scale associated with the
elastic deformation around the crack tip and the other associated with the
dissipation of energy at the crack tip. The morphology of the crack and the
nature of the crack growth depends much on the length of the crack in
comparison to these two length scales. A modification of the classical
Griffith’s theory has also been proposed [11] considering the discrete nature
of crack propagation which suggests: $\sigma_{c}\sim(2Yg/\pi(c+l))^{1/2}$,
where $c$ is a measure of the average increment of the length of the crack as
the external stress is increased adiabatically. A similar variation of nominal
stress with the crack length has also been observed experimentally [12, 13,
14, 15, 16] as well as numerically [17]. Recently, the balance of the internal
energies is explored in the mean-field limit [18] as well as with Langevin
dynamics [19] during a crack propagation through heterogeneous systems. In our
model, as we will see later, the energy balance criterion is taken care of
through the interplay between local stress profile and fluctuation of
strengths from point to point.

We present here the findings from the numerical study of cracking in a pre-
cracked one-dimensional fiber bundle model (FBM) [20]. FBM is perhaps the
simplest model which has been vastly used as a prototype model of fracture in
heterogeneous materials. The model is guided by threshold-activated dynamics
and shows different aspects of non-equilibrium statistical mechanics. The
dynamics and avalanches associated with crack propagation have been explored
earlier in the context of FBM [21, 22, 23]. In the present article, we differ
from previously explored studies by inserting a crack initially inside the
bundle and seeing how it affects the strength of the bundle, avalanche
statistics, and spatial correlation during the failure process. The reason
behind choosing one dimension is here the local stress concentration is most
prominent. As we go to a higher dimension, the model approaches the mean-field
limit and the presence of the pre-existing crack will be less effective. Our
main findings are that we do see the signature of two length scales ($\xi$ and
$\xi^{\ast}$), associated with the length of the pre-existing crack, in the
context of crack propagation in the bundle. We can relate the length scales
with the microscopic entities like number density and size distributions of
the cracks in the bundle. Further, as long as the length $l$ of the pre-
existing crack is smaller than $\xi$, the critical stress $\sigma_{c}$ of the
bundle remains independent of $l$. For $l>\xi$,
$\sigma_{c}\sim\sigma_{0}/(1+l/\xi)^{\alpha}$, where $\sigma_{0}$ is the
critical stress of the bundle in the absence of any pre-existing crack.
$\xi^{\ast}$, on the other hand, is an extreme limit of $l$ beyond which no
crack, other than the pre-existing one, develops within the bundle. We show
that the exponent $\alpha$ and the lengths $\xi$, $\xi^{\ast}$ depend on
$\beta$ and $\gamma$, the only two parameters which characterize the fiber
bundle model. $\beta$ is the strength of disorder (corresponding to the
strengths of the fibers) and is the spread of the fiber strength distribution.
$\gamma$ is a measure of the range up to which the stress of a broken fiber is
redistributed. The onset of localization has been studied earlier in the
context of fiber bundle model by tuning the disorder strength and without the
presence of any pre-existing crack [24, 25, 26]. In the present article, the
dependence of $\xi$, $\xi^{\ast}$ and $\alpha$ on $\beta$ and $\gamma$ and the
consequence of intrinsic crack resistance in crack propagation in the bundle
are discussed.

##  2 Description of Fiber Bundle Model

The fiber bundle model consists of a number of fibers (Hookean springs) held
perpendicularly between two parallel bars. The bars are pulled apart with
external stress $\sigma$ (force per fiber). The disorder is introduced in the
model as strengths of individual fibers which are randomly assigned following
a certain distribution. When the applied stress on a fiber crosses its
strength threshold, the fiber breaks irreversibly. The stress of the broken
fiber is then redistributed according to a certain stress redistribution
scheme. After redistribution, there might be further breaking of fibers due to
local enhancement of stress profile. The cascade of failure may continue
leading to global failure or stops if the stresses on fibers cannot reach
their threshold limit. The applied stress is to be increased to break the next
weakest fiber and the process is repeated until all fibers break.

We have considered a fiber bundle model where the strength of an individual
fiber is chosen randomly from a threshold distribution. We present here the
results for a power-law threshold distribution.

$p(\sigma_{\rm th})=\begin{cases}\displaystyle\frac{1}{2\beta\ln
10}\sigma_{\rm th}^{-1},&(10^{-\beta}\leq\sigma_{\rm th}\leq 10^{\beta})\\\
0&({\rm otherwise})\end{cases}$ (1)

Here $\beta$ denotes the width of the distribution or the strength of the
disorder. The motivation behind choosing such a distribution is the fact that
long-tailed distributions like Weibull [27] or power-law [28] has already been
observed for the distribution of material strength. We put a crack in a one-
dimensional FBM by removing $l$ consecutive fibers. If $a$ is the distance
between two consecutive fibers, the crack-length due to removal of $l$ fibers
will be $la$. We can set $a=1$ without losing any generality. From now on we
will address $l$ as the length of the crack that is inserted in the bundle.
This creates stress concentration on other unbroken fibers. The stress
concentration on the fibers will depend on the scheme one is adopting for
stress redistribution. We will mostly consider the local load sharing (LLS)
scheme (except for the results shown in figure 5) where the stress on the
broken fiber is redistributed on the nearest neighboring unbroken fibers.

Figure 1: The figures shows the stress profile of a 1d FBM of size $L$ under
local load sharing scheme (LLS) with a tensile force $F$ applied on it. A
crack of length $l(=3)$ (shown by the dotted red lines) is inserted in the
middle. This is equivalent to a crack of length $la$ (or $l$ if we consider
$a=1$). The local stress of all fibers are $\sigma(=F/L)$ except for the two
fibers at the crack tip that carries a stress $5\sigma/2$ each (see equation
2).

In this situation, if the bundle is subjected to an applied stress $\sigma$,
the stress on the two fibers at the crack tip due to a crack will be:

$\displaystyle\sigma_{t}(l)=\sigma\left(1+\displaystyle\frac{l}{2}\right)$ (2)

while a stress $\sigma$ is experienced by all other intact fibers. Figure 1
shows the local stress profile for a FBM of size $L=15$ with an external force
$F$ applied on it. This will create a stress $\sigma=F/L$ on the bundle in
absence of any pre-existing crack. A crack of size $l=3$ is inserted in the
middle by removing 3 consecutive fibers. This creates a local stress
$5\sigma/2$ at the crack tips while the stress on other fibers are $\sigma$.
If the stress on any fiber is higher than its threshold value, the fiber
breaks irreversibly and the stress of the broken fiber is redistributed
between the two nearest intact fibers at the two sides of it. For such
redistribution, we have adopted an algorithm that eliminates the memory effect
(similar to the algorithm in [29]). The algorithm we used is as follows. If a
fiber, carrying stress $\sigma_{b}$, is broken then after redistribution the
left and the right nearest neighboring unbroken fibers will experience the
stresses $\sigma_{l}$ and $\sigma_{r}$ respectively as

$\displaystyle\sigma_{l}\rightarrow\sigma_{l}+\displaystyle\frac{d_{r}}{d_{l}+d_{r}}\sigma_{b}$
$\displaystyle\sigma_{r}\rightarrow\sigma_{r}+\displaystyle\frac{d_{l}}{d_{l}+d_{r}}\sigma_{b}$
(3)

$d_{l}$ and $d_{r}$ are the distances of the left and right neighboring intact
fibers respectively from the broken fiber.

After the redistribution, either more fibers break due to increased local
stress profile or the bundle comes to a stable state where the redistributed
stress cannot overcome the threshold value for any fiber. The external stress
is then increased till a fiber breaks. We define $n(i)=\sigma(i)/\sigma_{e}$,
where $\sigma(i)$ and $\sigma_{e}$ are the local stress on $i$th fiber and
applied stress on the bundle respectively. For an external stress increment of
$\Delta\sigma_{e}$, the local stress increment on the $i$th fiber will be

$\displaystyle\Delta\sigma(i)=n(i)\Delta\sigma_{e}$ (4)

For fibers away from any broken fiber, $n(i)$ is 1. While if a fiber has an
adjacent cluster of broken fibers, $n(i)$ is greater than 1 depending on the
size of the cluster. The final value of $\sigma_{e}$ at which the global
failure occurs is known as the critical stress $\sigma_{c}$ or strength of the
bundle.

We have also studied a more general stress redistribution scheme where the
stress of a broken fiber $i$ is redistributed on the rest of the unbroken
fibers as:

$\displaystyle\sigma(j)\rightarrow\sigma(j)+\displaystyle\frac{r_{ij}^{-\gamma}}{Z}\sigma(i)$
(5)

where $\sigma(j)$ is the net stress on fiber $j$ at a distance $r_{ij}$ from
the broken fiber $i$. $Z$ is the normalizing factor and given by
$Z=\displaystyle\sum_{k}r_{ik}^{-\gamma}$, where $k$ runs over all intact
fibers. The initial local stress profile also depends on how far a certain
fiber is from the assigned crack. Eq.2 in this case will be modified as
follows:

$\displaystyle\sigma(j)=\sigma_{e}\left(1+\displaystyle\frac{R^{-\gamma}}{Z^{\prime}}l\right)$
(6)

Where $R$ is the distance of certain fiber from the crack. This means $R=1$
for the fibers at the notches, 2 for the next set of fibers, and so on.
$Z^{\prime}$ has the form $\displaystyle\sum R^{-\gamma}$ where the sum runs
over all intact fibers (which means $L-l$ fibers). In both equations 5 and 6,
$\gamma$ is a measure of the extent of the stress redistribution range. A high
$\gamma$ value corresponds to the LLS scheme whereas a low $\gamma$
corresponds to mean-field (MF) or equal load sharing (ELS) limit where the
stress of the broken fiber is equally distributed among all other intact
fibers in the bundle. A recent study [24] shows the existence of a critical
value $\gamma_{c}$ ($\approx 4/3$) for 1d FBM around which this change in
behavior from ELS to LLS takes place. The same was observed in 2d FBM as well
with a different $\gamma_{c}$ ($\approx 2$) [30]. In the present article, in
addition to the variation in $\beta$, we have also studied the crack
propagation for various values of $\gamma$.

##  3 Numerical Results

We present a numerical study of a pre-cracked fiber bundle model of size
$10^{5}$ in one dimension for $10^{4}$ realizations and for a wide range of
disorder strength $\beta$, stress relaxation range $\gamma$ and length $l$ of
the pre-existing crack. We address the question: ‘how does the size of the
pre-assigned crack relate to the rupture behavior of the bundle ?’

Figure 2: (a) Variation of $\rho$ with $B$ for various $l$ values and for
$\beta=1.0$. The red dots stand for maximum value $\rho_{m}$ of the patch
density. (b) Variation of $\rho_{m}$ with $l$ for $\beta=1.0$. For
$l>\xi^{\ast}$, $\rho_{m}$ falls to $1/L$ (pure nucleating). For $l<\xi$,
$\rho_{m}$ remains constant independent of l. (c) Variation of $P_{n}$ and
$P_{scn}$ with $l$ for $\beta=1.0$. Three distinct regions I, II and III are
seen. Region I: $P_{scn}$ is zero and $P_{n}$ is small indicating that the
effect of the pre-existing crack is not observed here. Region II: $P_{n}$ is
unity but $P_{scn}$ is low indicating that many cracks are observed here but
the pre-existing crack nucleates to create global failure. Region III: Both
$P_{n}$ and $P_{scn}$ have assumed the value unity indicating only the pre-
existing crack propagates here prior to global failure. (d) Average crack size
$s_{av}$ and maximum crack size $s_{max}$, just before global failure, are
plotted against $l$ for $\beta=1.0$. For $l<\xi$, both $s_{av}$ and $s_{max}$
are independent of $l$. Above $l>\xi$, $s_{max}\sim l$ suggesting coalescence
of the cracks to the pre-existing crack leading to fracture of the bundle.
Beyond $l>\xi^{\ast}$, $s_{av}=s_{max}\sim l$, suggesting nucleation of the
pre-existing crack only.

We first present the results for the LLS fiber bundle model. In figure 2(a),
we plot the density $\rho$ of the cracks (clusters of broken fibers) within
the bundle versus the fraction of broken bonds $B$ as the bundle ruptures with
time (time can be translated in this model in terms of redistributing steps).
$\rho$, at a certain time $t$, is given by the ratio of the number of cracks
at that time to the size $L$ of the system. $B$ is an increasing function of
time as it is given by a number of broken fibers at a certain time normalized
by the size $L$ of the system. $\rho$ has a value 0 and $1/L$ close to $B=0$
and $B=1$ respectively. In the former case, all fibers are intact (zero
patches) and in the latter case, all fibers are broken creating a single
patch. At an intermediate $B$, $\rho$ reaches a maximum value ($\rho_{m}$)
denoted by red dots in the figure. Beyond this point, crack clusters start to
coalesce and $\rho$ goes on decreasing till it reaches $1/L$ when there is a
single crack in the bundle.

Figure 2(b) shows the variation of $\rho_{m}$ with $l$. $\rho_{m}$ starts at a
constant value for low $l$ and decreases gradually to $\sim 1/L$ for high $l$
values. We see the signature of two length scales here. The first one is
$\xi$: for $l<\xi$, $\rho_{m}$ is independent of $l$. This indicates that for
low values of $l<\xi$, the pre-assigned crack has no effect on the fracturing
of the bundle and the bundle rupture by the random breaking of the fibers in
the bundle. The fracture here is dominated by the disorder and the stress
concentration has no effect on the fracture process. The second length scale
is $\xi^{\ast}$: for $l>\xi^{\ast}$, $\rho_{m}=1/L$, suggesting a pure
nucleating growth where a single crack, the pre-existing one, grows rupturing
the bundle. In this region, the cracking is dominated solely by stress
concentration. Naturally, $\xi,\xi^{\ast}$ depend on $\beta$ which will be
discussed later.

Figure 2(c) shows the variation of the probabilities, $P_{n}$ and $P_{scn}$
with $l$ for $\beta=1.0$. $P_{n}$ is the probability that at the final stage
of fracture, the instability of the bundle is triggered by the growth of the
crack which contains the pre-assigned crack. This is done by monitoring the
last fiber $j$ that breaks and causes instability and calculating the
probability that fiber $j$ is part of the enlarged pre-existing crack.
$P_{scn}$, on the other hand, is the probability for single crack nucleation,
where only one crack, the pre-existing one grows and ruptures the bundle. For
low $l$ in the region I ($l<\xi$), the low value of $P_{n}$ tells us that the
pre-assigned crack has no effect on the rupturing of the bundle. This is the
region where disorder dominates the fracturing process and sets the scale
$\xi$. In region III ($l>\xi^{\ast}$), $P_{scn}$ attains unity which means
that statistically, only the pre-assigned crack grows in the bundle rupturing
it. No other crack develops in the system. This is the extreme situation where
the fracture is dominated solely by stress concentration at the two ends of
the pre-assigned crack. The stress concentration sets the scale $\xi^{\ast}$.
In region II, both these two scales compete. In this region, many other cracks
develop in the system so that $P_{scn}$ is less than unity, but these cracks
coalesce with the pre-assigned crack finally rupturing the bundle. Region II
is thus dominated by both nucleation and coalescence.

Figure 2(d) shows how the average crack size $s_{av}$ and the maximum crack
size $s_{max}$ varies with $l$ for $\beta=1.0$. The diagonal dotted line is
the locus of $y\propto x$ while the vertical dotted lines show the positions
of $\xi$ and $\xi^{\ast}$ respectively. For $l<\xi$, both $s_{av}$ and
$s_{max}$ show no variation with $l$. $s_{max}$ is larger than $l$ and
$s_{av}$ is much smaller than $l$. This suggests that there are many cracks
developed in the bundle and the one that finally brings in instability and
ruptures the bundle has nothing to do with the pre-existing crack. For
$\xi^{\ast}>l>\xi$ , $s_{max}\sim l$ while $s_{av}$ remains much lower that
$l$. This suggests that though there are large number of cracks present in the
system, the largest one has size proportional to $l$. Finally, for
$l>\xi^{\ast}$, $s_{av}=s_{max}\sim l$ when only the pre-existing crack grows
in the bundle and no-other crack generate.

Figure 3: (a) The figure shows the variation of the critical stress
$\sigma_{c}$ with $l$ for LLS stress redistribution scheme and for
$\beta=1.0$. For $l<\xi$, $\sigma_{c}$ is independent of $l$. For $l>\xi$,
$\sigma_{c}\sim l^{-\alpha}$, where $\alpha\simeq 0.85$ is an exponent. (b)
Avalanche size distribution $P(\Delta)$ versus $\Delta$ for different lengths
$l$ of pre-existing cracks for LLS fiber bundle model at $\beta=1.0$. The
inset shows the the exponential decrease of $P(\Delta)$ with $\Delta$ for
small avalanche size $\Delta$ in log-normal scale. (c) The average size
$\langle\Delta_{f}\rangle$ of the large avalanches is plotted against $l$. (d)
The characteristic length $\Delta_{0}$ which defines the exponential fall of
the distribution is plotted against $l$.

Now that we have identified the two length scales associated with the crack
development in the fiber bundle model and established links between these
length scales and the microscopic entities like size and number distribution
of clusters of broken fibers, we would like to discuss how the macroscopic
observables (often measured in cracking experiments) behave at different
length scales. To this end, we have measured the critical stress and the
avalanche behavior of the bundle for various $l$ values. These two quantities
have been discussed a lot in the context of crack propagation in materials
[31] and in fiber bundle model [20].

Figure 3(a) shows the variation of critical stress $\sigma_{c}$ with $l$ in
LLS model and for $\beta=1.0$. We find for $l<\xi$, $\sigma_{c}$ is
independent of $l$. For $l>\xi$, $\sigma_{c}$ decreases with $l$ as

$\displaystyle\sigma_{c}=\displaystyle\frac{\sigma_{0}}{\left(1+\displaystyle\frac{l}{\xi}\right)^{\alpha}}$
(7)

Here $\xi$ appears as the cut-off length and $\alpha$ is an exponent whose
value depends on $\beta$ and $\gamma$. We will discuss this dependence later.
As is mentioned in the introduction, this form of dependence of $\sigma_{c}$
on $l$ is very much seen in experiments [12, 13, 14, 15, 16] and simulations
[17] and also has been suggested theoretically [1, 2, 3]. It is to be noted
that even for $l=0$, there is a finite value of $\sigma_{c}$ unlike that in
Griffith’s law. This is a consequence of intrinsic crack resistance.

Figure 4: (a) The figure shows the variation of $\xi$ and $\xi^{\ast}$ with
$\beta$ for LLS FBM. We observe three regions: I, II and III, depending on the
spatial correlation during failure. Region I: conventional LLS limit where the
crack does not play any role. Region II: multiple cracks but only the pre-
existing crack nucleates and creates global failure. Region III: only a single
crack, the pre-existing one, propagates prior to global failure. No other
crack originates within the bundle. (b) Variation of $\alpha$ with $\beta$ for
LLS FBM. $\alpha$ has a high value, close to 1, for low $\beta$ and then
decreases with $\beta$. For $\beta>0.5$, $\alpha$ remains constant around 0.85
independent of the disorder. strength.

Figure 3(b) shows the avalanche size distribution for different lengths $l$ of
the pre-assigned crack in LLS fiber bundle at $\beta=1.0$. An avalanche is
defined as the number of fibers broken in between two consecutive stress
increments. The avalanches in the fiber bundle model show a robust behavior in
its size distribution [32, 33] as well as have information about upcoming
catastrophic events [34]. We have presented the results in terms of $\Delta$,
which is the avalanche size or the number of broken fibers in one avalanche
normalized by size $L$ of the system. Figure 3(b) shows the variation of
$P(\Delta)$ with $\Delta$. In abrupt brittle failure at low disorder (both ELS
[35] and LLS [36] FBM), there is no need for stress increment after the
weakest fiber is broken. In this limit, the bundle breaks in a single
avalanche. On the other hand, at a moderate disorder, the avalanche size
distribution shows a scale-free behavior with a universal exponent -5/2 [32]
and an exponential nature [33] for ELS and LLS scheme respectively. We observe
the same exponential decay of the avalanche size distribution
$P(\Delta)\sim\exp(-\Delta/\Delta_{0})$ for finite-size avalanches and
distribution for largest avalanche $\Delta_{f}$ prior to rupture even if there
is a pre-existing crack within the bundle. In the regime, $l<\xi$, the
characteristic length $\Delta_{0}$ for the exponential decay and the average
value $\langle\Delta_{f}\rangle$ of the largest avalanche remain almost
constant. However, both these quantities increase rapidly with $l$ for $l>\xi$
(see figure 3(c) & (d)). For large $l$, the rapid exponential fall of finite
avalanches and a large value of the large avalanche indicates that the bundle
ruptures by destabilization of the crack by large avalanches.

Figure 5: (a) Variation of $\sigma_{c}$ with $l$ for $\gamma=0.5$
($<\gamma_{c}$) and $\gamma=4.0$ ($>\gamma_{c}$) and for $\beta=1.0$. The
behavior matches with Eq.7 (LLS FBM) and Eq.11 (MF FBM) when $\gamma$ is 4.0
and 0.5 respectively. (b), (c) and (d) respectively shows the variation of
$\xi$, $\xi^{\ast}$ and $\alpha$ with $\gamma$. High $\gamma$ corresponds to
the LLS limit. On the other hand, when $\gamma$ is low, both $\xi$ and
$\xi^{\ast}$ increases with decreasing $\gamma$. This is due to the fact that
in low $\gamma$, the local stress concentration does not play a crucial role
and the pre-existing crack will be visible to the model only when the length
of the crack is large. $\alpha$ decreases with decreasing $\gamma$ suggesting
less effect of the crack on the critical stress.

We would now discuss how the length scales and the exponent $\alpha$ depend on
$\beta$. Figure 4(a) shows the variation of above mentioned two characteristic
lengths, $\xi$ and $\xi^{\ast}$, when the strength of disorder $\beta$ is
increased. We observe three distinct regions. Region I: $l<\xi$ and the pre-
assigned crack has no effect on the cracking of the bundle. The critical
stress is independent (see figure 3) of the crack length. The dynamics of the
system in absence of any pre-assigned crack has been explored earlier in
detail [37]. This region is dominated by the disorder and stress
redistribution has not much effect on cracking. For high $\beta$, the rupture
sequence of the fibers is random, whereas, for low $\beta$, cracks nucleate at
suitable soft spots independent of the pre-assigned crack. Region II: In this
region, disorder competes with stress redistribution in cracking. Failure is
dominated by crack length. The critical stress falls with $l$ in a scale-free
manner with an exponent $\alpha$ (see figure 3). The probability that the pre-
existing crack causes final instability in the bundle is unity (see figure
2c), which means the pre-existing crack nucleates though there are other
cracks developed within the bundle. Region III: This region is dominated by
stress redistribution. This is an extreme limit of the crack-length dominated
region. Here only the pre-existing crack propagates within the bundle and no
other crack is originated. The critical stress follows the scale-free decay
with $l$ here as well. One should notice that for low $\beta$, $\xi$ and
$\xi^{\ast}$ coincide with each other. In this limit, the disorder strength is
not strong enough to trap the pre-existing crack unless the length of the
crack is small itself. As a result, we either see random failure events
(region I) at low $l$ or pure nucleating failure beyond a certain length scale
(region III). The increase in length scales with $\beta$ makes sense as a
large crack has to be inserted within the bundle at high disorder so that the
stress concentration overcomes the high fluctuation in threshold strength
values, in other words, the intrinsic crack resistance.

Figure 4(b) shows the variation of $\alpha$ with $\beta$ for LLS FBM. $\alpha$
has a high value 1 for $\beta=0$, which is consistent with equation 2 since
all fibers have the same threshold and the entire bundle ruptures in a single
avalanche as soon as the weakest fiber is broken. As $\beta$ increases,
$\alpha$ starts to decrease and remains constant around 0.85 independent of
the disorder beyond $\beta=0.5$.

So far, we have presented our findings for a particular value of $\beta$ and
in the LLS fiber bundle model. In principle, we have two parameters in the
model: the disorder measured by its strength $\beta$ and the range of stress
redistribution characterized by $\gamma$. We have seen that the origin of two
length scales is essentially due to the competition between stress
concentration on the fibers due to stress redistribution in the cracking
process and the disorder in the threshold strengths of the fibers. We would
expect that our results will also depend on the stress redistribution range
characterized by $\gamma$ in addition to the disorder strength characterized
by the width of the threshold distribution of the fibers.

We now present our results for the variation of $\xi$, $\xi^{\ast}$ and
$\alpha$ with $\gamma$. We have already discussed that a low value of $\gamma$
corresponds to the MF or ELS version of the model. The critical stress, in
presence of a pre-existing crack, can be calculated in the MF limit
analytically. We will then compare the analytical result with the numerical
behavior of critical stress with $l$ keeping $\gamma=0.5$ (which is much below
$\gamma_{c}=4/3$ for 1d FBM [24]). For fiber strength distribution as given by
(1), the cumulative distribution is given by

$\displaystyle
P(\sigma_{th})=\displaystyle\int_{10^{-\beta}}^{\sigma_{th}}p(\sigma_{th})d\sigma_{th}=\displaystyle\frac{\ln\sigma_{th}}{2\beta\ln
10}+\displaystyle\frac{1}{2}$ (8)

where $p(\sigma_{th})=\sigma_{th}^{-1}/(2\beta\ln 10)$ is the normalized
density function. For applied force $F$ and an extension $\epsilon$ of the
fibers, functional form between $F$ and $\epsilon$ can be expressed as

$\displaystyle
F(\epsilon)=L\epsilon\left(\displaystyle\frac{1}{2}+\displaystyle\frac{\ln\epsilon}{2\beta\ln
10}\right)$ (9)

The function $F(\epsilon)$ has a maximum at critical extension
$\epsilon_{c}=10^{\beta}/e$ ($e\approx 2.72$), satisfying
$dF(\epsilon)/d\epsilon=0$. The critical value $\sigma_{t}(l)$ required to
break the model in presence of the crack of length $l$ is generated by a
critical external stress $\sigma_{c}$,

$\displaystyle\sigma_{t}(l)=\sigma_{c}\left(1+\displaystyle\frac{l}{L-l}\right)=\displaystyle\frac{F_{c}}{L}=\displaystyle\frac{F(\epsilon_{c})}{L}=\displaystyle\frac{10^{\beta}}{2\beta
e\ln 10}$ (10)

The critical stress $\sigma_{c}$ is given as follows

$\displaystyle\sigma_{c}=\displaystyle\frac{10^{\beta}}{2\beta e\ln
10}\left(1-\displaystyle\frac{l}{L}\right)$ (11)

This shows that the critical stress increases with disorder strength $\beta$
and decreases with crack-length $l$.

Figure 5 shows how the cracking behavior depends on the stress release range
$\gamma$. Figure 5(a) shows for a high value of $\gamma$ (= 4.0), $\sigma_{c}$
shows the same behavior with $l$ as in the LLS fiber bundle model:
$\sigma_{c}$ is independent of $l$ for $l<\xi$ and decays with a power law for
$l>\xi$. The behavior for $\sigma_{c}$ at high $\gamma$ is fitted with Eq.7
with $\xi=28$ and $\alpha=0.85$. Both $\xi$ and $\alpha$ are functions of
$\gamma$. For a low $\gamma$ (= 0.5) value, the model is in the MF (or ELS)
limit and we observe $\sigma_{c}$ versus $l$ behavior to be very close to the
mean-field expression given by Eq.11.

Figure 5(b), (c) and (d) show the variations of $\xi$, $\xi^{\ast}$ and
$\alpha$ respectively with $\gamma$. In the limit of high $\gamma$, all the
three parameters $\xi$, $\xi^{\ast}$ and $\alpha$ saturates towards the values
we obtained for the LLS fiber bundle model. As $\gamma$ decreases both $\xi$
and $\xi^{\ast}$ increases. This happens as with decreasing $\gamma$ the
effect of stress concentration at the notches of the pre-existing crack
becomes less effective and a larger crack is required to make the stress
concentration effective again. This in turn increases the sizes of the
characteristic length scales which finally diverges around $\gamma_{c}$ (the
dotted lines in figure 5(b) and (c)). $\alpha$ on the other hand decreases as
$\gamma$ decreases up to $\gamma\geq\gamma_{c}$. This is also due to the fact
that the stress concentration becomes less effective as $\gamma$ increases and
that increases the chances of more resistance to cracking. This lowers the
value of $\alpha$. Due to the same reason, earlier $\alpha$ was observed to
decrease with an increase in $\beta$ as well. Below $\gamma_{c}$, as the model
is in the mean-field limit, the behavior of $\sigma_{c}$ is not scale-free and
$\alpha$ does not have meaning there.

Figure 6: (a) $\sigma_{c}$ vs $l$ for system sizes ranging in between $10^{3}$
and $10^{5}$. For $l<\xi$, $\sigma_{c}$ is independent of $l$ but decreases
with increasing $L$ (see inset (i)). The system size effect is not visible in
the region $l>\xi$ (see inset (ii)). (b) Variation of characteristic length
$\xi$ with disorder strength $\beta$ for system sizes ranging in between
$10^{3}$ and $10^{5}$. $\xi$ is not observed to change with size of the
system. We set $\beta=1.0$.

We briefly discuss any possible effect of system size on our results. Figure
6(a) shows the variation of $\sigma_{c}$ with $l$ for system sizes $10^{3}$,
$5\times 10^{3}$, $10^{4}$, $5\times 10^{4}$ and $10^{5}$ with LLS scheme.
$\beta$ is kept constant at 1.0. When $l<\xi$, a reduction in critical stress
is observed with the increasing size of the bundle (see the inset (i) of the
same figure). This reduction in $\sigma_{c}$ is observed due to the weakest
link of chain like failure dynamics in this model as the weakest link itself
scales to lower values as the size of the bundle increases [20]. On the other
hand, we do not see any system size dependence $\xi$ for $l>\xi$ (see the
inset (ii) of the same figure). The probable reason for this is in this limit
the critical stress is already lower than the value we observe due to the
weakest link effect. This happens due to the fact that a very small external
stress can create large stress at the crack notches. The redistributed stress
depends on the crack-length and disorder strength and not on the system size.
Though we don’t claim the same size-effect of $\sigma_{c}$ at higher
dimensions since the system size scaling rule of critical stress is not valid
there [29]. Figure 6(b) explicitly shows that the size of the system does not
have any effect on this characteristic length $\xi$. Both the regions,
disorder-dominated and crack-dominated, remain unchanged when we increase the
system size and we can conclude both these regions will survive in the
thermodynamic limit as well.

The universal behavior of our results has been explored briefly for uniform
and Weibull distribution for fiber strengths. The preliminary study does not
show any dependence of the results on the nature of the distribution. However,
this aspect needs to be studied in detail.

##  4 Discussion

In conclusion, we have mimicked the strength experiments generally done by
stressing a specimen with a pre-existing crack by a numerical study of
stressing fiber bundle model with a pre-existing crack in it. The important
observations from cracking in engineering samples are the pronounced effect of
lattice trapping or crack resistance, modification of Griffith’s law, and
occurrence of two length scales arising from the competition between the large
scale elastic stress relaxation and short scale energy dissipation near the
crack tip. The fiber bundle model is perhaps the simplest model for fracture
in heterogeneous materials having essentially two parameters: the disorder in
the fiber strength thresholds and regime over which the stress of a broken
fiber is redistributed. We clearly see the existence of two length scales
associated with the cracking of the bundle. We relate the two length scales
with the microscopic observables like the number and size distribution of the
cracks in the bundle. We see the modification of Griffith’s law similar to
what is found in engineering specimens and we determine the range of validity
of the law in the case of the fiber bundle model. The other spectrum of this
is the high disorder scenario where the failure events are random in space.
This has been observed numerically in 2d random register network [38, 39] as
well as in random spring network [40]. Finally, we have determined the
dependence of these two length scales and the exponent that characterizes the
modified Griffith’s law on the strength of the disorder and the range of
stress redistribution from a failed fiber in the bundle.

##  5 Acknowledgements

SR was supported by the Research Council of Norway through its Centres of
Excellence funding scheme, project number 262644. TH was supported by Japan
Society for the Promotion of Science (JSPS) Grants-in-Aid for Scientific
Research (KAKENHI), Grants Nos. 16H06478 and 19H01811. PR acknowledges JSPS
Invitational Fellowships for Research in Japan and Earthquake Research
Institute, Tokyo for infrastructural support for the visit.

## References

  * [1] Lawn, B. R., & Wilshaw T. R. Fracture of Brittle Solids. Cambridge University Press, Cambridge (1975).
  * [2] Knott, J. F. Fundamentals of Fracture Mechanics. Butterworths (1973).
  * [3] Broberg, K. B. Cracks and Fracture. Academic Press, New York (1999).
  * [4] N. Bernsteina, and D. W. Hess, Phys Rev Lett 91, 025501 (2003).
  * [5] A. Mattoni, L. Colombo, and F. Cleri, Phys Rev Lett 95, 115501 (2005).
  * [6] W. A. Curtin, J Mater Res 5, 1549 (2000).
  * [7] R. Perez and P. Gumbsch, Phys Rev Lett 84, 5347 (2000).
  * [8] J. R. Rice, J Mech Phys Solids 40, 239 (1992).
  * [9] R. Thomson, C. Hsieh, and V. Rana, J Appl Phys 42, 3154 (1971).
  * [10] R. Long, C. Y. Hui, J. P. Gong, and E. Bouchbinder, Ann. Rev. of Cond. Matt. Phys. 12, 71 (2021).
  * [11] M. P. Wnuk, and A. Yavari, Eng. Fract. Mech. 75, 1127 (2008).
  * [12] Z. P. Bazant, Metal Journal of Eng. Mec., Vol. 110, 518-535 (1984).
  * [13] Z. P. Bazant and J. Planas, Fracture and size effect in concrete and other quasibrittle materials. CRCPress (1997).
  * [14] A. Carpinteri, Theo. and App. Frac. Mech. 2 39 (1984).
  * [15] K. Duan, X. Z. Hu and F. H. Wittmann, Materials and Structures/Materiaux et Constructions 36, 74 (2003).
  * [16] R. W. Armstrong, Phil. Trans. R. Soc. A 373, 20140124 (2014).
  * [17] T. Nojima, and H. Tsuyoshi, J. Soc. Mat. Sci., Japan 44(499), 451-456 (1995).
  * [18] H. B. da Rochaa, and L. Truskinovsky, J. Mec. Phys. of Sol. 158, 104646 (2022).
  * [19] Y. Feng, J. Fan, and E. B. Tadmor, J. Mec. Phys. of Sol. 159, 104715 (2022).
  * [20] A. Hansen, P. C. Hemmer, and S. Pradhan, The Fibre Bundle Model. Weinheim, Germany: Wiley-VCH (2015).
  * [21] A. Delaplace, S. Roux, G. Pijaudier-Cabot, Journal of engineering mechanics, 127(7), 646 (2001).
  * [22] S. Patinet, D. Vandembroucq, and A. Hansen, Eur. Phys. J. Spec. Top. 223, 2339 (2014).
  * [23] F. Villette, J. Baroth, F. Dufour, J. F. Bloch, and S. R. Du Roscoat, International Journal of Fracture, 221(1), 87 (2020).
  * [24] S. Biswas, S. Roy, and P. Ray, Phys Rev E 91, 050105(R) (2015).
  * [25] S. Roy, Front. Phys. 9, 752086 (2021).
  * [26] A. Stormo, K. S. Gjerden, and A. Hansen, Phys. Rev. E 86, 025101(R) (2012).
  * [27] D. Munz and T. Fett, Ceramics: Mechanical Properties, Failure Behavior, Materials Selection. Springer-Verlag Berlin Heidelberg (1999).
  * [28] R. A. Fisher, and L. H. C. Tippett, Proc. Cambridge Philos. Soc. 24: 180–190 (1928).
  * [29] S. Sinha, J. T. Kjellstadli and A. Hansen, Phys Rev E 92, 020401(R) (2015).
  * [30] R. C. Hidalgo, Y. Moreno, F. Kun, and H. J. Herrmann, Phys. Rev. E 65, 046148 (2002).
  * [31] S. Biswas, P. Ray, and B. K. Chakrabarti, Statistical Physics of Fracture Breakdown and Earthquake, Wiley-VCH, Berlin Germany (2015).
  * [32] P. C. Hammer and A. Hansen, ASME J Appl Mech 59, 909 (1992).
  * [33] M. Kloster, A. Hansen, and P. C. Hemmer, Phys Rev E 56, 2615 (1997).
  * [34] V. Kadar, G. Pal, and F. Kun, Scientific Reports 10, 2508 (2020).
  * [35] S. Roy, and P. Ray, Europhysics Letters 112, Number 2, 26004 (2015).
  * [36] S. Roy, Phys Rev E 96, 042142 (2017).
  * [37] S. Roy, S. Biswas, and P. Ray, Phys Rev E 96, 063003 (2017).
  * [38] A. A. Moreira, C. L. N. Oliveira, A. Hansen, N. A. M. Araujo, H. J. Herrmann, and J. S. Andrade Jr., Phys. Rev. Lett. 109, 255701 (2012).
  * [39] A. Shekhawat, S. Zapperi, and J. P. Sethna, Phys. Rev. Lett. 110, 185505 (2013).
  * [40] P. Ray and G. Date, Physica A: Statistical Mechanics and its Applications 229, 26 (1996).


# Traversable wormholes in $f(R)$ gravity sourced by a cloud of strings

Parangam Goswami ID E-mail<EMAIL_ADDRESS>Department of Physics,
Assam University, Cachar - 788011, Assam, India  Anshuman Baruah ID E-mail:
<EMAIL_ADDRESS>Department of Physics, Assam University, Cachar -
788011, Assam, India  Atri Deshamukhya ID Correspondence to:
<EMAIL_ADDRESS>Department of Physics, Assam University, Cachar -
788011, Assam, India

###### Abstract

Wormhole solutions in General Relativity (GR) require exotic matter sources
that violate the null energy condition (NEC), and it is well known that
higher–order modifications of GR and some alternative matter sources can
support wormholes. In this study, we explore the possibility of formulating
traversable wormholes in $f(R)$ modified gravity, which is perhaps the most
widely discussed modification of GR, with two approaches. First, to
investigate the effects of geometrical constraints on the global
characteristics, we gauge the $rr$–component of the metric tensor, and employ
Padè approximation to check whether a well–constrained shape function can be
formulated in this manner. We then derive the field equations with a
background of string cloud, and numerically analyse the energy conditions,
stability, and amount of exotic matter in this space–time. Next, as an
alternative source in a simple $f(R)$ gravity model, we use the background
cloud of strings to estimate the wormhole shape function, and analyse the
relevant properties of the space–time. These results are then compared with
those of wormholes threaded by normal matter in the simple $f(R)$ gravity
model considered. The results demonstrate that wormholes with NEC violations
are feasible; however, the wormhole space-times in the simple $f(R)$ gravity
model are unstable.

Keywords— Wormhole, string cloud, modified gravity, energy conditions

###### Contents

  1. 1 Introduction
  2. 2 Traversable wormholes
    1. 2.1 A novel shape function
  3. 3 Traversable wormholes in $f(R)$ gravity
    1. 3.1 Cloud of Strings as a source
    2. 3.2 Traversable wormholes in $f(R)=\alpha R^{m}-\beta R^{n}$ gravity with a string cloud background
    3. 3.3 Traversable wormholes in $f(R)=\alpha R^{m}-\beta R^{n}$ gravity with ordinary matter
  4. 4 Discussion and Conclusion

##  1 Introduction

The governing equation of the theory of General Relativity (GR) is the
Einstein field equation (EFE), and solutions to this set of coupled
differential equations have been remarkably successful in accounting for
phenomena ranging from black holes [1] to the evolution of the universe.
Spherical symmetry is an important physically relevant constraint in solving
the EFEs, and wormholes are an intriguing prospect among spherically symmetric
solutions, explored since the early days of GR [2, 3]. Ellis [4] and Bronnikov
[5] independently reported the first traversable Lorentzian wormhole solution,
and the geometric requirements were established in detail by Morris & Thorne
in 1988 [6]. Wormholes are exact solutions to the EFEs, and can be interpreted
as bridges connecting two different asymptotically flat regions of space–time
via a _throat_. The topology of the wormhole interior is non-trivial, while
the topology at the boundaries remains simple [7]. While theoretical as of
now, wormholes are of crucial importance in fundamental physics, especially
considering quantum entanglement [8] and quantum gravity. Moreover, it has
recently been shown that wormholes may actually mimic astrophysical black
holes in some observations [57]. A fundamental problem in such space–times is
that for the _throat_ to remain open for signal propagation, the behavior of
the matter sources supporting such a space–time is _exotic_ in that the null
energy condition (NEC) is violated [6]. Specifically, this implies that
observers in an inertial frame measure negative energy densities at the
throat, which is unphysical. Such behaviour can be avoided in modified gravity
theories.

Despite its success, GR fails to explain cosmological phenomena such as
late–time accelerated cosmic expansion [9, 10] and the so called ‘inflationary
epoch’ [11]. Modified theories of gravity address these shortcomings [12, 13,
14], and the extra degrees of freedom in modified gravity theories also
enables one to evade NEC violations in wormhole space–times. Moreover,
wormhole solutions are an inherent feature of most modifications of GR. One of
the simplest modifications to GR is $f(R)$ modified gravity, where the Ricci
scalar $R$ in the Einstein–Hilbert action is replaced by some arbitrary
function of it [15]. This simple modification to GR can lead to a host of
models that can independently meet both cosmological and solar system tests
[16]. It has been shown that the extra degrees of freedom in $f(R)$ gravity
arising from higher order curvature terms may lead to scenarios where the
matter content satisfies the NEC in wormhole space–times [17, 18]. Wormholes
in the framework of $f(R)$ gravity have been studied extensively in different
iterations of $f(R)$ gravity (for instance, in Refs. [19, 20, 21, 22, 23, 24,
25, 26, 27, 28]).

Another inherent limitation of GR is that it provides only a classical
description of gravity, and the theory is non–renormalisable at high energy
(small length) scales. String theory [29] is perhaps the strongest contender
to a unified paradigm of gravity with the other fundamental forces, and it
posits that the fundamental constituent of matter and energy are extended
objects, instead of point–like ones. Precisely, the extended objects are
considered as one–dimensional relativistic strings, and the interactions of
strings on a classical level provide better models of several fundamental
interactions [30, 31, 32]. To this end, string clouds as a gravitational
source have been studied extensively in literature. A general solution to the
EFEs for a spherically symmetric cloud of strings was first reported in [33],
with an emphasis on the energy conditions. Properties of compact objects such
as black holes in the background of a string cloud have been reported
previously [34, 35, 36]. Traversable wormholes in the background of a string
cloud have been reported [37], with an emphasis on the amount of exotic matter
and stability of the wormhole configuration against radial perturbations.
Recently, a detailed study of the properties of traversable wormholes
surrounded by a string cloud in the framework of $f(R)$ gravity have been
reported with analysis of the quasi normal modes (QNMs) of the wormhole
solutions [38].

In this study, we investigate traversable wormholes in $f(R)$ gravity with two
motivations: First, to investigate the effects of geometrical constraints, we
gauge the $rr$–component of the metric tensor, and employ Padé approximation
to check whether a well–constrained shape function can be formulated in this
manner. We then derive the field equations with the background of a string
cloud, and numerically analyse the energy conditions, stability, and amount of
exotic matter in this space–time. Next, we use a background of string cloud to
estimate the wormhole shape function in a simple $f(R)$ gravity model, with
$f(R)=\alpha R^{m}-\beta R^{n}$ [39], and analyse the relevant properties of
the space–time as in the previous case. These results are then compared with
those of wormholes threaded by normal matter in the same modified gravity
model considered.

The remainder of this manuscript is organised as follows. In Sec. 2, we
discuss the traversable wormhole geometry. A novel shape function is proposed
using Padé approximation in Sec. 2.1. In Sec. 3, we present the modified EFEs
in the general framework of $f(R)$ gravity. In Sec. 3.1 we analyse the various
energy conditions, the stability in terms of TOV equation, and the amount of
exotic matter required to sustain traversible wormholes in the framework of
$f(R)$ gravity with a background of string cloud. Next, in Sec. 3.2 we
estimate the shape function in the considered form of $f(R)=\alpha R^{m}-\beta
R^{n}$ gravity model [39] with a string cloud background and analyse the
properties of the wormhole space-time. For a comparative analysis, in Sec. 3.3
we present wormhole solutions supported by ordinary matter for the $f(R)$
model [39] considered in Sec. 3.2 using our proposed shape function obtained
by employing Padé approximation in Sec. 2.1. Finally, we conclude the work
with some remarks in Sec. 4. We adhere to use the natural system of units
($G=c=1$) throughout the work.

##  2 Traversable wormholes

Morris & Thorne [6] used the following metric ansatz to describe a static,
spherically symmetric space-time

$ds^{2}=-e^{2\Phi(r)}dt^{2}+\frac{dr^{2}}{1-\frac{b(r)}{r}}+r^{2}{d{\theta}}^{2}+r^{2}sin^{2}{\theta}d{\phi}^{2}$
(1)

Eq. (1) is the line element of a traversable wormhole. The proper radial
coordinate $l(r)=\int_{r_{0}}^{r}\frac{dr}{\sqrt{1-\frac{b}{r}}}$ should be
well behaved throughout the space-time in order to avoid singularities. It
imposes the constraint $\frac{b}{r}\leq 1$ at the throat. This space-time is a
special case of the Ellis–Bronnikov [4, 5] space-time described in terms of
$l(r)$ as

$ds^{2}=-dt^{2}e^{2\Phi(l)}+dl^{2}+r^{2}(l)({d{\theta}}^{2}+r^{2}sin^{2}{\theta}d{\phi}^{2})$
(2)

where the throat is located at some minimum $r(l)$.
The metric function $\Phi(r)$ in Eq. (1) is known as the red-shift function,
and the first co-efficient of the line element in Eq. (1) provides a measure
of the gravitational red-shift. The topological configuration of the space-
time is determined by the second coefficient of the line element in Eq. (1),
and the metric function $b(r)$ is known as the shape function. At some minimum
value of the radial co-ordinate $r$, the throat of the wormhole is located at
some arbitrary value $r_{0}$. A significant aspect of traversable wormholes is
that the throat should not be surrounded by an event horizon. Horizons in
spherically symmetric space-times are described by physically non-singular
surfaces at $g_{00}=-e^{2\Phi}\rightarrow 0$, and this results in the
constraint that throughout the space-time $\Phi(r)$ should be well defined.
Moreover, the geometric constraints on the shape function $b(r)$ demanded by
traversability are: (i) $b(r_{o})=r_{o}$, (ii)
$\frac{b(r)-b^{\prime}(r)r}{b^{2}}>0$, (iii) $b^{\prime}(r_{o})-1\leq 0$, (iv)
$\frac{b(r)}{r}<1,\forall r>r_{o}$, (v) $\frac{b(r)}{r}\rightarrow 0$ as
$r\rightarrow\infty$, where prime denotes a derivative with respect to the
radial co-ordinate $r$. The energy density $\rho$, radial pressure $p_{r}$,
and transverse pressure $p_{t}$ of the matter sources are constrained by these
conditions on the metric functions through the EFEs. Therefore, while
constructing traversable wormhole configurations violations of the energy
conditions appear owing to these constraints. Next, we propose a novel form
the shape function and check the viability of the shape function by analysing
the various constraints. In addition, in this work we consider tideless
traversable wormhole solutions described by a constant red-shift function
$\Phi^{\prime}(r)=0$.

###  2.1 A novel shape function

We examine the following functional form as a probable shape function

$\displaystyle
b(r)=r_{0}\left[\log\left(\frac{r}{r_{0}}\right)+\coth(r_{0})\tanh(r)\right]^{a}$
(3)

where, $r_{0}$ is the throat, and $a$ is a free parameter. The viability of
the shape function $b(r)$ depends on several constraints as discussed in Sec
2. Analysing these constraints for the functional form of Eq. (3), insights
can be obtained about the throat radius $r_{0}$, and the parameter $a$.
However, it can be observed Eq. (3), does not satisfy the condition of
$b(r_{0})=r_{0}$, which is an important consequence of traversable wormholes,
the other required conditions for the viability of the shape function, such as
the asymptotic flatness, and the flaring out condition are also not satisfied.
Thus, in order to have a plausible form of $b(r)$, the next step involves
adopting the Padé approximation for the functional form in Eq. (3), and to
check if a viable form of $b(r)$ can be obtained that satisfies all the
required constraints on a shape function. Use of rational expansions made by
Padé functions is common in existing literature [40, 41, 42, 43]. The Padé
approximation is built on the Taylor series expansion. For a given function
such as $f(z)=\sum_{i=1}^{\infty}c_{i}z^{i}$, expanding the series with the
coefficients $c_{i}$, the $(n,m)$ Padé approximant ratio is given as [44],

$\displaystyle
P_{n,m}(z)=\frac{\sum_{k=0}^{n}a_{k}z^{k}}{1+\sum_{\sigma=0}^{m}b_{\sigma}z^{\sigma}}$
(4)

The most simple forms of Padé approximant are of the orders $(1,0)$ & $(0,1)$.
Thus, considering the $(1,0)$ approximation and expanding Eq. (3) about the
throat radius $r_{0}$, the shape function is

$\displaystyle b(r)=r_{0}-a(r_{0}-r)\left[1+r_{0}\leavevmode\nobreak\
\text{cosech}(r_{0})\text{sech}(r_{0})\right]$ (5)

With the form as described in Eq. (5), the various conditions required for the
viability are analysed and it becomes evident that the function in Eq. 5
satisfies all the necessary conditions to be a shape function for $r_{0}=3$
and $0<a<1$. Moreover, the function satisfies the condition $b(r_{0})=r_{0}$.
We therefore propose $b(r)=r_{0}-a(r_{0}-r)\left[1+r_{0}\leavevmode\nobreak\
\text{cosech}(r_{0})\text{sech}(r_{0})\right]$ as a new viable shape function.

The plots of the asymptotic flatness $b(r)/r\rightarrow 0$ as
$r\rightarrow\infty$, and flaring out condition
$\frac{b(r)-b^{\prime}(r)r}{b^{2}}>0$ are shown in Figure 1.

Figure 1:  Profile of the (a) Asymptotic flatness $\frac{b(r)}{r}\rightarrow
0$ as $r\rightarrow\infty$, and (b) Flaring out condition
$\frac{b(r)-b^{\prime}(r)r}{b^{2}}>0$ vs. $r$ for $0<a<1$

With the obtained viable shape function, next the traversable wormhole
solutions are analysed in the framework of $f(R)$ gravity with a background of
string cloud.

##  3 Traversable wormholes in $f(R)$ gravity

A general form of the action in $f(R)$ modified theories is [19]

$S=\int d^{4}x\sqrt{-g}\left[f(R)+\mathcal{L}_{m}\right]$ (6)

Using the metric formalism of $f(R)$ gravity, the modified EFEs are obtained
as

$FR_{\mu\nu}-\frac{1}{2}f(R)\,g_{\mu\nu}-\nabla_{\mu}\nabla_{\nu}F+g_{\mu\nu}\Box
F=\,T^{m}_{\mu\nu}\,,$ (7)

where $F\equiv df/dR$, and
$T^{m}_{\mu\nu}=\frac{-2}{\sqrt{-g}}\frac{\delta\mathcal{L}_{m}}{\delta
g^{\mu\nu}}$ is the stress-energy tensor of background matter. We consider
that the spherically symmetric space-time is represented by the line-element
in Eq. (1). Contracting Eq. (7) yields:

$FR-2f(R)+3\Box F=T$ (8)

Here, $T$ is the trace of the stress-energy tensor of matter and $\Box F$ is:

$\displaystyle\Box
F=\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}g^{\mu\nu}\partial_{\nu}F)=\left(1-\frac{b}{r}\right)\left[F^{\prime\prime}-\frac{b^{\prime}r-b}{2r^{2}(1-b/r)}\,F^{\prime}+\frac{2F^{\prime}}{r}\right]$
(9)

with $F^{\prime}=df(R)/dR$ and $b^{\prime}=d\,b(r)/dr$. Substituting Eq. (8)
in Eq. (7) yields the modified EFEs as:

$G_{\mu\nu}\equiv R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=T^{{\rm eff}}_{\mu\nu}\,$
(10)

Here, $T^{{\rm eff}}_{\mu\nu}$ is the effective stress-energy tensor,
responsible for the energy condition violations and is generally interpreted
as a gravitational fluid. $T^{{\rm eff}}_{\mu\nu}$ comprises of the matter
stress energy tensor $T^{m}_{\mu\nu}$ and the curvature stress-energy tensor
$T^{{\rm c}}_{\mu\nu}$ given by

$\displaystyle
T^{c}_{\mu\nu}=\frac{1}{F}\left[\nabla_{\mu}\nabla_{\nu}F-\frac{1}{4}g_{\mu\nu}\left(RF+\Box
F+T\right)\right]$ (11)

Assuming that the geometry of the wormhole is threaded by an anisotropic
distribution of matter

$T_{\mu\nu}=(\rho+p_{t})U_{\mu}\,U_{\nu}+p_{t}\,g_{\mu\nu}+(p_{r}-p_{t})\chi_{\mu}\chi_{\nu}\,,$
(12)

where $U^{\mu}$ is the four-velocity, and $\chi^{\mu}$ represents a unit
space-like vector.

With the line element in Eq. (1), the modified EFEs can be expressed as the
following [19]

$\displaystyle\rho=\frac{Fb^{\prime}}{r^{2}}$ (13) $\displaystyle
p_{r}=-\frac{bF}{r^{3}}+\frac{F^{\prime}}{2r^{2}}(b^{\prime}r-b)-F^{\prime\prime}\left(1-\frac{b}{r}\right)$
(14) $\displaystyle
p_{t}=-\frac{F^{\prime}}{r}\left(1-\frac{b}{r}\right)+\frac{F}{2r^{3}}(b-b^{\prime}r)$
(15)

The Ricci scalar is given as $R=\frac{2b^{\prime}}{r^{2}}$. With the explicit
forms of energy density $\rho$, radial pressure $p_{r}$, and transverse
pressure $p_{t}$ in Eqs. (13)-(15), the various energy conditions, stability
and the amount of exotic matter required for the wormhole configuration can be
analysed.

###  3.1 Cloud of Strings as a source

A cloud of strings is analogous to the perfect fluid models of gas and dust.
However, the difference is that it comprises one-dimensional objects extended
along some specific direction. The string cloud can exist in different
geometrical configurations such as planar, axisymmetric, or spherical [45]. A
general solution for the spherical distribution of the cloud of strings was
reported in [45]. In addition, thermodynamic properties of string gas has been
reported in [46]. For a spherically symmetric string cloud in four-dimensions,
the energy momentum tensor has the following non-null components [45]

${T^{t}}_{t}={T^{r}}_{r}=-\frac{\eta^{2}}{r^{2}}$ (16)

where $\eta$ is a constant related to the total energy of the string cloud.
With the form of the shape function as described in Eq. (5), and using Eq.
(13), we get

$F(r)=-\frac{\eta^{2}}{a\left(r_{0}\leavevmode\nobreak\
\text{cosech}(r_{0})\leavevmode\nobreak\ \text{sech}\leavevmode\nobreak\
(r_{0})+1\right)}$ (17)

With the obtained $F(r)$ and the field equations Eqs. (13)-(15), the numerical
analyses are conducted to obtain the energy conditions, check the stability,
and estimate the amount of exotic matter.

For any observer traversing a time-like curve to detect the energy density of
the matter field to be positive, the stress-energy tensor of matter must
adhere to some sets of inequalities known as the energy conditions in GR [47,
18]. The weak energy condition (WEC) implies $\rho\geq 0$, the NEC implies
$\rho+p_{r}\geq 0$, and $\rho+p_{t}\geq 0$, whereas the strong energy
condition (SEC) implies $\rho+p_{r}+2p_{t}\geq 0$. The profile of energy
conditions with a background of string cloud is presented with the throat
radius fixed at $r_{0}=3$. Emphasis is made on the effect of the string cloud
constant $\eta$ on the energy conditions. Results are discussed by fixing the
parameter $a=0.5$. For the whole range of $a$, $0<a<1$, the energy conditions
show similar behaviour. Previous studies indicate that the constant associated
with the total energy of the string cloud should be small [33, 37].

Figure 2 shows the profile of the NEC terms $\rho+p_{r}$ and $\rho+p_{t}$. It
can be observed that the first NEC term is satisfied at the throat $r_{0}=3$
for the considered values of the $\eta$ as
$\eta=0.1,\,0.3,\,\text{and}\leavevmode\nobreak\ 0.5$. However, the second NEC
term $\rho+p_{t}$ is violated for all the values of $\eta$. Owing to the
violation of the second NEC term, as a whole the NEC is inferred to be
violated.

Figure 2:  Profile of the NEC terms (a) $\rho+p_{r}$, and (b) $\rho+p_{t}$
respectively vs. $r$ with $b(r)$ as in Eq. 5

From figure 3, it can be observed that the WEC is violated at the wormhole
throat for all the values of $\eta$. It is also evident from the fact that the
non-null components of the stress-energy tensor of the string cloud comes with
a minus sign as shown in Eq. (16). Moreover, it is seen that the SEC exhibits
an oscillatory (indeterminate) behaviour along the radial co-ordinate $r$.
With increasing $r$, the oscillation between positive and negative values
decreases; however, the behaviour extends asymptotically.

Figure 3:  Profile of the (a) WEC $\rho$, and (b) SEC $\rho+p_{r}+2p_{t}$ vs.
$r$ respectively with $b(r)$ as in Eq. 5

In addition to the energy conditions, analysing two other parameters, viz. the
equation of state (EoS) parameter $\omega=p_{r}/\rho$, and anisotropy
parameter $\Delta=p_{t}-p_{r}$ turns out to be useful. Information about the
nature of the matter source threading the wormhole geometry can be obtained
from the EoS parameter, and the attractive or repulsive nature of the space-
time geometry (geometrical viability) can be understood by the anisotropy
parameter.

Figure 4:  Profile of the (a) EoS parameter $\omega$, and (b) anisotropy
parameter $\Delta$ vs. $r$ respectively with $b(r)$ as in Eq. 5

Figure 4 shows the behaviour of the EoS and the anisotropy parameters. Near
the wormhole throat, the EoS parameter is $\omega<-1$, signifying a phantom-
like behaviour of the string cloud for all values of $\eta$. The anisotropy
parameter $\Delta<0$ near the wormhole throat, signifying an attractive nature
of the space-time geometry. The summary of the energy conditions are presented
in Table 1.

Table 1: Summary of the energy conditions discussed in Sec. 3.1 Terms |  Result | Interpretation
---|---|---
$\rho$ |  $<0$ near throat, for $\eta=0.1,0.3,0.5$ |  WEC violated at throat
$\rho+p_{r}$ |  $>0$ near throat, for $\eta=0.1,0.3,0.5$ |  NEC satisfied at throat
$\rho+p_{t}$ |  $<0$ near throat, for $\eta=0.1,0.3,0.5$ |  NEC violated at throat
$\rho+p_{r}+2p_{t}$ |  oscillates, for $\eta=0.1,0.3,0.5$ |  SEC indeterminate
$\omega$ |  $<-1$ near throat, for $\eta=0.1,0.3,0.5$ |  phantom-like source at throat
$\Delta$ |  $<0$ near throat, for $\eta=0.1,0.3,0.5$ |  attractive geometry at throat

First reported in the context of neutron stars [48, 49], the Tolman-
Oppenheimer-Volkov (TOV) equation provides information regarding the stability
of stellar structures. In order to probe the stability of wormholes in terms
of the hydrostatic, gravitational, and anisotropic forces in the space-time, a
more generalized version of the formalism was developed in [50]. The
generalized TOV equation [50, 51] is given as

$-\frac{dp_{r}}{dr}-\frac{\epsilon^{\prime}(r)}{2}(\rho+p_{r})+\frac{2}{r}(p_{t}-p_{r})=0,$
(18)

where $\epsilon(r)=2\Phi(r)$. $F_{\mathrm{h}}$ represents the hydrostatic
force, $F_{\mathrm{g}}$ the gravitational force, and $F_{\mathrm{a}}$, the
anisotropic force. These three terms of the TOV equation can determine the
equilibrium anisotropic mass distribution [51] in that stable stellar
structures satisfy Eq. (18).

$F_{\mathrm{h}}=-\frac{dp_{r}}{dr},\;\;\;\;\;\;\;\;F_{\mathrm{a}}=\frac{2}{r}(p_{t}-p_{r}),\;\;\;\;\;\;\;\;F_{\mathrm{g}}=-\frac{\epsilon^{{}^{\prime}}}{2}(\rho+p_{r})$
(19)

Owing to the constant red-shift function $\Phi^{\prime}(r)=0$, the
gravitational force is $F_{\mathrm{g}}=0$ in the analysis.

Using the averaged null energy condition,
$\int_{\lambda_{1}}^{\lambda_{2}}T_{ij}k^{i}k^{j}d\lambda\geq 0$, evaluated
along the radial coordinate $r$, the amount of exotic matter in wormhole
space-times can be estimated. However, the amount of energy condition
violating matter can be estimated in a more generalised manner by using a
volume integral instead of the line integral, namely the volume integral
quantifier (VIQ) [52, 53, 54], which is defined as

$I_{v}=\oint[\rho+p_{r}]dV=8\pi\int_{r_{0}}^{s}(\rho+p_{r})r^{2}dr$ (20)

After matching the wormhole space-time with an exterior metric and cutting off
the stress-energy tensor at some $r=s$ away from the throat, an estimate of
the amount of NEC violating matter can be obtained by using the VIQ. For
arbitrarily small quantities of NEC violating matter, the requirement is
$I_{v}\rightarrow 0$ as $s\rightarrow r_{0}$ [52, 53].

Figure 5: Profile of the (a) $F_{\mathrm{h}}$, and $F_{\mathrm{a}}$ vs. $r$
(b) VIQ with $b(r)$ as in Eq. 5

Figure 5 shows the terms of the TOV equation and the profile of the VIQ. It is
seen that that the hydrostatic force $F_{\mathrm{h}}$, and the anisotropic
force $F_{\mathrm{a}}$ cancel each other asymptotically, signifying a stable
configuration. In addition, it can be observed that the VIQ, $I_{v}\rightarrow
0$ as $s\rightarrow r_{0}$, indicates that the wormhole solution is feasible
with arbitrarily small amounts of exotic matter. The terms of the TOV and VIQ
equations are shown with $\eta=0.5$. For $\eta=0.1,\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ 0.3$ the corresponding terms of the TOV and VIQ
equations depict a similar behaviour.

Thus, wormhole solutions can be formulated in a model independent manner using
the novel shape function, Eq. (5), in the framework of $f(R)$ gravity with a
string cloud background. The next motive remains to check whether a feasible
form of the shape function can be obtained in a particular form of $f(R)$
gravity model with the string cloud cloud background.

###  3.2 Traversable wormholes in $f(R)=\alpha R^{m}-\beta R^{n}$ gravity with
a string cloud background

In this section, we obtain a form of the shape function $b(r)$ by fixing
$f(R)$ and analyse its viability. For the analysis, we consider a simple form
of $f(R)$ gravity model as [39]

$f(R)=\alpha R^{m}-\beta R^{n}$ (21)

where $\alpha$ and $\beta$ are positive constants and $m$ and $n$ are positive
integers satisfying the condition $m>n$ [39]. Using Eq. (16) in Eq. (13), we
obtain

$-\frac{\eta^{2}}{r^{2}}=\frac{Fb^{\prime}}{r^{2}}$ (22)

Considering the $f(R)$ model as described in Eq. (21), and obtaining $F$, the
shape function can be found from Eq. (22). For simplifying the calculations,
the integers $m$ and $n$ are set as $m=2$ and $n=1$ satisfying the condition
$m>n$. Again, using the expression for the Ricci scalar
$R=\frac{2b^{\prime}}{r^{2}}$, $F$ is found in terms of $r$. Eq. (22) reduces
to a quadratic equation in $b^{\prime}$ of the form

$4\alpha{b^{\prime}}^{2}-b^{\prime}\beta r^{2}+\eta^{2}r^{2}=0$ (23)

Solving Eq. (23) and considering only the positive root for mathematical
feasibility, we get

$b^{\prime}=\frac{\beta
r^{2}+\sqrt{\beta^{2}r^{4}-16\alpha\eta^{2}r^{2}}}{8\alpha}$ (24)

Integrating Eq. (24), we get

$b=\frac{1}{8\alpha}\left[\frac{\beta
r^{3}}{3}+\frac{\left(\beta^{2}r^{2}-16\alpha\eta^{2}\right)\sqrt{\beta^{2}r^{4}-16\alpha\eta^{2}r^{2}}}{3\beta^{2}r}\right]+c$
(25)

where $c$ is a constant of integration. In order to evaluate the constant $c$,
we use the condition at the wormhole throat, $b(r_{0})=r_{0}$. This leads to
the following form of the shape function

$\displaystyle b$
$\displaystyle=r_{0}+\frac{1}{24\alpha\beta^{2}}\left[\frac{\beta^{3}r^{4}+\left(\beta^{2}r^{2}-16\alpha\eta^{2}\right)\sqrt{\beta^{2}r^{4}-16\alpha\eta^{2}r^{2}}}{r}\right.$
$\displaystyle-\left.\frac{\beta^{3}r_{0}^{4}+\left(\beta^{2}r_{0}^{2}-16\alpha\eta^{2}\right)\sqrt{\beta^{2}r_{0}^{4}-16\alpha\eta^{2}r_{0}^{2}}}{r_{0}}\right]$
(26)

With the shape function in Eq. (3.2), the various conditions required for the
viability of the shape function as described in Sec 2 are analysed. It is
interesting to note that the asymptotic flatness $\frac{b}{r}\rightarrow 0$ as
$r\rightarrow\infty$, and the flaring out condition
$\frac{b-b^{\prime}r}{b^{2}}>0$ are satisfied only for negative values of the
constants $\alpha$ and $\beta$ ($\alpha=-0.5$ and $\beta=-0.8$). The profile
of the asymptotic flatness and flaring out condition are shown in Figure 6
with different values of $\eta$.

Figure 6: Profile of the (a) asymptotic flatness $\frac{b}{r}$ and (b) flaring
out condition $\frac{b-b^{\prime}r}{b^{2}}$ vs. $r$ with $r_{0}=1$,
$\alpha=-0.5$ and $\beta=-0.8$, for $b(r)$ in Eq. 3.2

With the shape function in Eq. (3.2), and the $f(R)$ model described by Eq.
(21), the energy conditions, stability and the amount of exotic matter
required for the wormhole configuration are analysed next.

Using Eqs. (13)-(15), the various energy conditions are analysed with the
throat at $r_{0}=1$, and by fixing the model parameters of the $f(R)$ gravity
model (Eq. (21)) as, $m=2,\,n=1$, $\alpha=-0.5$ and $\beta=-0.8$.

Figure 7: Profile of the NEC terms (a) $\rho+p_{r}$ and (b) $\rho+p_{t}$ vs.
$r$ with $b(r)$ as in Eq. 3.2 Figure 8: Profile of the (a) WEC $\rho$ and (b)
SEC $\rho+p_{r}+2p_{t}$ vs. $r$ with $b(r)$ as in Eq. 3.2 Figure 9: Profile of
the (a) EoS parameter $\omega$ and (b) anisotropy parameter $\Delta$ vs. $r$
with $b(r)$ as in Eq. 3.2

Figure 7 shows the NEC terms $\rho+p_{r}$, and $\rho+p_{t}$. It can be
observed that the first NEC term $\rho+p_{r}$ is violated at the wormhole
throat. The second NEC term $\rho+p_{t}$ is satisfied at the wormhole throat.
However, owing to violation to the first NEC term, as a whole the NEC is
inferred to be violated.

Figure 8 shows the WEC and the SEC. It can be seen that the WEC is violated at
the wormhole throat, and this violation can again be attributed to the
negative sign associated with the non-null components of the stress-energy
tensor of the string cloud as in Eq. (16). The SEC is satisfied at the throat
and also asymptotically, which is again an interesting point to note as in
$f(R)$ gravity, the SEC should be asymptotically violated to account for the
late-time accelerated expansion of the universe.

Table 2: Summary of the energy conditions discussed in Sec. 3.2 Terms |  Result |  Interpretation
---|---|---
$\rho$ |  $<0$ near throat, for $\eta=0.1,0.3,0.5$ |  WEC violated at throat
$\rho+p_{r}$ |  $<0$ near throat, for $\eta=0.1,0.3,0.5$ |  NEC violated at throat
$\rho+p_{t}$ |  $>0$ near throat, for $\eta=0.1,0.3,0.5$ |  NEC satisfied at throat
$\rho+p_{r}+2p_{t}$ |  $>0$ near throat, for $\eta=0.1,0.3,0.5$ |  SEC satisfied at throat
$\omega$ |  $>0$ near throat, for $\eta=0.1,0.3,0.5$ |  normal matter like source at throat
$\Delta$ |  $>0$ near throat, for $\eta=0.1,0.3,0.5$ |  repulsive geometry at throat
Figure 10: Profile of the (a) $F_{\mathrm{h}}$, and $F_{\mathrm{a}}$ vs. $r$
and (b) VIQ with $b(r)$ as in Eq. 3.2

Figure 9 shows the variation of the EoS parameter $\omega$ and the anisotropy
parameter $\Delta$. It can be seen that the EoS parameter $\omega$ is
$\omega>0$ near the wormhole throat for all values of $\eta$, signifying that
the string cloud as the source behaves like ordinary matter without any
phantom-like behaviour. The anisotropy parameter $\Delta$ is positive near the
wormhole throat for all the values of $\eta$, signifying a repulsive geometry
at the throat. The summary of the energy conditions are presented in Table 2.

Figure 10 shows the terms of the TOV equation and the VIQ. It is seen that the
hydrostatic force $F_{\mathrm{h}}$, and anisotropic force $F_{\mathrm{a}}$
does not cancel each other out, signifying that the wormhole configuration is
unstable. From the VIQ, it is evident that $I_{v}\rightarrow 0$ as
$s\rightarrow r_{0}$, indicating that the wormhole configuration can be
obtained with arbitrarily small amounts of exotic matter.

With the results presented in Sec 3.1 and Sec 3.2, it is clear that wormhole
solutions in $f(R)$ gravity with background a string cloud is feasible with
characteristic violation of the NEC. However, with $f(R)=\alpha R^{m}-\beta
R^{n}$ gravity, the wormhole configuration is not stable. In order to have a
better understanding, the next section presents results of wormhole solution
in $f(R)=\alpha R^{m}-\beta R^{n}$ gravity supported by ordinary matter and
with the form of the shape function as in Eq. (5).

###  3.3 Traversable wormholes in $f(R)=\alpha R^{m}-\beta R^{n}$ gravity with
ordinary matter

With the shape function in Eq. (5), the field equations Eqs. (13)-(15) are
analysed with the form of the $f(R)$ model as in Eq. (21). We assume that the
wormhole geometry is threaded by ordinary matter with the stress-energy tensor
$T^{\mu}_{\nu}=\text{daig}[-\rho(r),p_{r}(r),p_{t}(r),p_{t}(r)]$. Owing to the
properties of the shape function, the throat radius is fixed at $r_{0}=3$, and
the free parameter $a$ is considered for the whole range $0<a<1$. The model
parameters of the $f(R)=\alpha R^{m}-\beta R^{n}$ gravity are fixed as,
$\alpha=0.8$, $\beta=0.5$, and $m=2,\,n=1$. The various energy conditions,
stability of the wormhole space-time and the VIQ are presented below.

Figure 11 shows the profile of the NEC terms $\rho+p_{r}$ and $\rho+p_{t}$. It
is evident that that the first NEC term $\rho+p_{r}$ is violated at the
wormhole throat. However, the second NEC term $\rho+p_{t}$ is satisfied at the
throat. Owing to the violation of the first NEC term, as a whole the NEC is
considered to be violated.

Figure 12 depicts the WEC and the SEC. It is seen that the WEC is violated
near the wormhole throat. The SEC is satisfied at the wormhole throat and
asymptotically for the whole range of $a$, $0<a<1$. It is interesting to note
that even with ordinary matter as the source, the SEC is satisfied
asymptotically for $f(R)=\alpha R^{m}-\beta R^{n}$ gravity.

Figure 11: Profile of the NEC terms (a) $\rho+p_{r}$ and (b) $\rho+p_{t}$ vs.
$r$ with $b(r)$ as in Eq. 5 Figure 12: Profile of the (a) WEC $\rho$ and (b)
SEC $\rho+p_{r}+2p_{t}$ vs. $r$ with $b(r)$ as in Eq. 5

Figure 13 shows the variation of the EoS parameter and the anisotropy
parameter. It can be observed that as expected, the EoS parameter is
$\omega>0$ near the wormhole throat signifying that the source threading the
wormhole geometry has no phantom-like behaviour. The anisotropy parameter
$\Delta>0$ near the throat, signifying a repulsive geometry at the throat. The
energy conditions are summarised in Table LABEL:tab3.

Figure 13: Profile of the (a) EoS parameter $\omega$ and (b) anisotropy
parameter $\Delta$ vs. $r$ with $b(r)$ as in Eq. 5

Figure 14 shows the corresponding terms of the TOV equation, and the VIQ. It
is seen that both the hydrostatic force $F_{\mathrm{h}}$, and the anisotropic
force $F_{\mathrm{a}}$ does not cancel each other out asymptotically for the
whole range of $0<a<1$, signifying that the wormhole configuration is not
stable. From the VIQ it is evident that $I_{v}\rightarrow 0$ as $s\rightarrow
r_{0}$, indicating that the wormhole solution can be obtained with arbitrarily
small amount of exotic matter. However, the VIQ is only shown for $a=0.5$, and
for the whole range of $0<a<1$, the VIQ has a similar profile.

Figure 14: Profile of the (a) $F_{\mathrm{h}}$, and $F_{\mathrm{a}}$ vs. $r$
and (b) VIQ with $b(r)$ as in Eq. 5

##  4 Discussion and Conclusion

In this study, we investigated the traversable wormholes in the framework of
$f(R)$ gravity with a background of string clouds. A novel form of the shape
function using Padé approximation was proposed for the analysis. Using this
shape function a specific form of $F(r)$ was obtained in Eq. (17) in order to
analyse the EFEs. However, it is observed that the $F(r)$ depends on the
parameters of the metric functions and on the string cloud parameter $\eta$.
Eliminating these dependence and obtaining a cosmologically viable form of
$f(R)$ remains an open issue. The results demonstrate that stable wormhole
solutions with characteristic violation of the energy conditions (especially
the NEC) are feasible in the framework of $f(R)$ gravity with the new shape
function in the background of a string cloud. In addition, the EoS parameter
near the wormhole throat indicates that the string cloud has phantom-like
properties. The anisotropy parameter $\Delta<0$ near the throat signifying an
attractive geometry. In addition, the SEC shows an indeterminate behaviour as
it oscillates between positive and negative values. Therefore, to have a
better understanding of the wormhole solutions, traversable wormholes are
considered in a string cloud background with a simple form of $f(R)$ gravity
model, $f(R)=\alpha R^{m}-\beta R^{n}$. It is interesting to note that the
shape function in this case can yield viable wormhole geometries for negative
values of $\alpha$ and $\beta$. The results demonstrate characteristic
violations of the NEC. Interestingly, the SEC is satisfied both at the throat
and asymptotically. The EoS parameter $\omega>0$ at the throat and indicates
that that the string cloud behaves as normal matter. The anisotropy parameter
$\Delta>0$ near the wormhole throat and signifies a repulsive geometry at the
throat. However, the TOV equation shows that the wormhole space-time is
unstable. For a comparative analysis, wormholes threaded by normal matter are
analysed in $f(R)=\alpha R^{m}-\beta R^{n}$ gravity, using the Padé
approximate shape function given in Eq. (5). The constants $\alpha$ and
$\beta$ are considered with positive values. The results demonstrate
characteristic violation of the NEC and shows that the wormhole space-time is
unstable. The SEC is satisfied at the wormhole throat and asymptotically.

Although wormholes have not been detected yet, studying these exact solutions
of the EFEs provide with far reaching insight into the nature of space-time
and fundamental building blocks of the universe. Wormholes also play a crucial
role in several important issues such as the cosmic censorship conjecture
[55]. Studies have reported strong indications for their presence in the form
of black hole mimickers [56, 57]. Further studies regarding the observational
constraints of these wormhole solutions remain an open issue to be addressed
in the near future.

## References

  * [1] K. Schwarzschild, Preuss. Akad. Wiss. Berlin (Math. Phys.) 1916, 189–196 (1916)
  * [2] L. Flamm, Phys. Z. 17, 448 (1916)
  * [3] A. Einstein, N. Rosen, Phys. Rev. 48(1), 73 (1935)
  * [4] H.G. Ellis, J. Math. Phys. 14(1), 104–118 (1973)
  * [5] K.A. Bronnikov, Acta Phys. Pol p. B4 (1973)
  * [6] M.S. Morris, K.S. Thorne, Am. J. Phys. 56, 395–412 (1988)
  * [7] M. Visser, _Lorentzian wormholes: From Einstein to Hawking_ (American Institute of Physics Melville, NY, USA, 1996)
  * [8] L. Susskind, Fortschritte der Phys. 64(1), 72–83 (2016)
  * [9] S. Perlmutter, M.S. Turner, M. White, Phys. Rev. Lett. 83(4), 670 (1999)
  * [10] A.G. Riess, P.E. Nugent, R.L. Gilliland, B.P. Schmidt, J. Tonry, M. Dickinson, R.I. Thompson, T. Budavari, S. Casertano, A.S. Evans, et al., Astrophys. J. 560(1), 49 (2001)
  * [11] A. Linde, _Particle physics and inflationary cosmology_ , vol. 5 (CRC press, London, 1990)
  * [12] A.A. Starobinsky, Phys. Lett. B 91(1), 99–102 (1980)
  * [13] E.J. Copeland, M. Sami, S. Tsujikawa, Int. J. Mod. Phys. D 15(11), 1753–1935 (2006)
  * [14] S. Capozziello, M. De Laurentis, Phys. Rep. 509(4-5), 167–321 (2011)
  * [15] T.P. Sotiriou, V. Faraoni, Rev. Mod. Phys. 82(1), 451 (2010)
  * [16] J.Q. Guo, Solar system tests of $f(R)$ gravity.  Int. J. Mod. Phys. D 23(04), 1450,036 (2014)
  * [17] S. Capozziello, F.S. Lobo, J.P. Mimoso, Phys. Rev. D 91(12), 124,019 (2015)
  * [18] A. Baruah, A. Deshamukhya, J. Phys. Conf. Ser. 1330(1), 012,001 (2019)
  * [19] F.S. Lobo, M.A. Oliveira, Phys. Rev. D 80(10), 104,012 (2009)
  * [20] N. Furey, A. DeBenedictis, Class. Quantum Gravity 22(2), 313 (2004)
  * [21] N. Godani, G.C. Samanta, Mod. Phys. Lett. A 34(28), 1950,226 (2019)
  * [22] N. Godani, G.C. Samanta, Int. J. Mod. Phys. D 28(02), 1950,039 (2019)
  * [23] T. Azizi, Int. J. Theor. Phys. 52, 3486–3493 (2013)
  * [24] A.K. Mishra, U.K. Sharma, V.C. Dubey, A. Pradhan, Astrophys. Space Sci. 365(2), 34 (2020)
  * [25] C.G. Boehmer, T. Harko, F.S. Lobo, Phys. Rev. D 85(4), 044,033 (2012)
  * [26] G. Mustafa, Z. Hassan, P. Sahoo, Ann. Phys. 437, 168,751 (2022)
  * [27] A. Baruah, P. Goswami, A. Deshamukhya, Int. J. Mod. Phys. D p. 2250119 (2022)
  * [28] A. Baruah, P. Goswami, A. Deshamukhya, New Astron. 99, 101,956 (2023)
  * [29] S. Mukhi, Class. Quantum Gravity 28(15), 153,001 (2011)
  * [30] M. Kalb, P. Ramond, Phys. Rev. D 9(8), 2273 (1974)
  * [31] P.S. Letelier, Phys. Rev. D 15(4), 1055 (1977)
  * [32] F. Lund, T. Regge, Phys. Rev. D 14(6), 1524 (1976)
  * [33] P.S. Letelier, Phys. Rev. D 20(6), 1294 (1979)
  * [34] S.G. Ghosh, S.D. Maharaj, Phys. Rev. D 89(8), 084,027 (2014)
  * [35] A. Belhaj, Y. Sekhmani, Gen. Relativ. Gravit. 54(2), 17 (2022)
  * [36] D.V. Singh, S.G. Ghosh, S.D. Maharaj, Phys. Dark Universe 30, 100,730 (2020)
  * [37] M.G. Richarte, C. Simeone, Int. J. Mod. Phys. D 17(08), 1179–1196 (2008)
  * [38] D.J. Gogoi, U.D. Goswami, J. Cosmol. Astropart. Phys. 2023(02), 027 (2023)
  * [39] S. Nojiri, S.D. Odintsov, Phys. Rev. D 77(2), 026,007 (2008)
  * [40] S. Capozziello, R. D’Agostino, O. Luongo, Int. J. Mod. Phys. D 28(10), 1930,016 (2019)
  * [41] C. Gruber, O. Luongo, Phys. Rev. D 89(10), 103,506 (2014)
  * [42] Y.N. Zhou, D.Z. Liu, X.B. Zou, H. Wei, Eur. Phys. J. C 76, 1–13 (2016)
  * [43] S. Capozziello, O. Luongo, L. Mauro, Eur. Phys. J. Plus 136, 1–14 (2021)
  * [44] H. Padé, in _Ann. Sci. de l’Ecole Norm. Superieure_ , vol. 9 (1892), pp. 3–93
  * [45] J.M. Graça, I.P. Lobo, I.G. Salako, Chin. Phys. C 42(6), 063,105 (2018)
  * [46] T.G. Mertens, H. Verschelde, V.I. Zakharov, J. High Energy Phys 2014(11), 1–35 (2014)
  * [47] E. Curiel, in _Towards a theory of spacetime theories_ (Springer, New York, NY, 2017), pp. 43–104
  * [48] R.C. Tolman, _Relativity, thermodynamics, and cosmology_ (The Clarendon Press, Oxford, London, 1934)
  * [49] R.J. Oppenheimer, G.M. Volkoff, Phys. Rev. 55(4), 374 (1939)
  * [50] V. Gorini, U. Moschella, A.Y. Kamenshchik, V. Pasquier, A.A. Starobinsky, Phys. Rev. D 78(6), 064,064 (2008)
  * [51] J.P. de Leon, Gen. Relativ. Gravit. 25(11), 1123–1137 (1993)
  * [52] M. Visser, S. Kar, N. Dadhich, Phys. Rev. Lett. 90(20), 201,102 (2003)
  * [53] S. Kar, N. Dadhich, M. Visser, Pramana 63(4), 859–864 (2004)
  * [54] F.S.N. Lobo, F. Parsaei, N. Riazi, Phys. Rev. D 87(8), 084,030 (2013)
  * [55] R. Penrose, Gen. Relativ. Gravit. 34(7), 1141–1165 (2002)
  * [56] R.N. Izmailov, A. Bhattacharya, E.R. Zhdanov, A.A. Potapov, K.K. Nandi, Eur. Phys. J. Plus 134(8), 384 (2019)
  * [57] K.K. Nandi, R.N. Izmailov, A.A. Yanbekov, A.A. Shayakhmetov, Phys. Rev. D 95(10), 104,011 (2017)