1 Introduction

A large number of papers on static traversable wormholes have been written in the last decades [2, 11, 16, 20, 24, 26, 31, 52, 53]. In fact, wormholes are bridges between two branes, universes, or just connections of two points at the manifold. Generally wormholes are asymptotically flat. There were many proposals for the wormhole models. One of the first models was proposed by [15], and it was called the Einstein–Rosen bridge. Einstein–Rosen bridges are vacuum solutions of Einstein field equations, and this type of wormholes are just an internal part of the maximally extended Schwarzschild black hole metric. A maximally extended metric means that this metric has no boundaries and the geodesic lines of the particle can be laid infinitely far into the future. So if the spacetime is maximally extended, then there must be present the so-called white hole interior. The exterior of the white hole is often called another universe. The white hole and the second universe are needed in order to extend the trajectory of a particle that fell beyond the event horizon of Schwardschild’s black hole infinitely far into the future.

In Figure. 1 one can see the Penrose diagram, where \(i^0\) is the infinitely far spacelike point, \(i^-\) is the infinitely distant past, \(i^+\) is the infinitely distant future point. Thus timelike curves lie from \(i^-\) to \(i^+\). Therefore, similarly to the light cone, here \(\mathscr {I}^+\) and \(\mathscr {I}^-\) are lightlike infinitely distant future/past. The upper shaded part of the figure is the interior of our universe black hole, with a singularity at \(r=0\), and the bottom shaded triangle is respectively the white hole interior of another universe with a singularity at \(r=0\). As well, \({\mathscr {H}}^+\) is the black hole horizon and \({\mathscr {H}}^-\) is the white hole antihorizon, \(\Sigma \) is a spacelike geodesic trajectory through both universes (Cauchy surface). Here both universes are just Minkowski manifolds.

After the theoretical prediction of Einstein–Rosen bridges and Schwarzschild wormholes, many astrophysicists and cosmologists have begun to search for the possibility of the existence of traversable wormholes. One of the first wormhole options, and at the moment one of the most plausible, is the option proposed by [31]. This is the static traversable Morris-Throne wormhole. This type of wormhole can connect two points of spacetime, and its throat is located in the bulk (with \(d_{\mathrm {bulk}}>4\)). This exact solution is a good definition for a traversable wormhole, but as it turned out, this solution implies the presence of an exotic matter at the throat. Therefore, the Null Energy Condition (NEC) is violated in classical GR gravity. By varying parameter values in modified theories of gravity, we could solve this problem, or at least minimize the amount of NEC violating matter.

1.1 Chosen modified theories of gravity to research

Although the general theory of relativity was and is a wonderful theory that describes our universe well enough, it still has its problems. General theory of relativity (further – GR) is a non-renormalizable theory of gravity, and therefore, can not be conventionally quantized [48]. This is just one of many problems of classical gravity in general relativity. Moreover, GR fails to explain the recent cosmological observations [14]. In order to overpass these problems, astrophysicists introduced in the literature modified theories of gravity. The novelty of the modified theories of gravity is that new geometrodynamic terms are introduced in the gravitational field by the modification of the Einstein–Hilbert (EH) Action Integral. In the following, we briefly introduce the modified theories of gravity of our consideration.

Fig. 1
figure 1

Penrose diagram for maximally extended Schwarzschild black hole metric (code provided by Robert McNees)

1.2 \(f({\mathscr {R}})\) gravity

\(f({\mathscr {R}})\) gravity is the typical and most popular choice of modified gravity theory, which modifies Einstein–Hilbert action and replaces Ricci scalar in the EH action with arbitrary function of Ricci scalar \(f({\mathscr {R}})\). The theory was originally proposed in [6]. It has drawn the attention of cosmologists because it can provide a geometric mechanism for the description of inflation [5, 21, 49] and of the dark energy problem [7, 33].

1.3 Metric-Palatini gravity

Metric-Palatini is a completely different type of modified gravity theory. In Metric-Palatini gravity in addition to Ricci scalar in EH action function \(f(\mathscr {\hat{R}})\) is introduced, which is an arbitrary function for the Palatini scalar, which is constructed from the metric tensor and the Levi-Civita connection. Hybrid Metric-Palatini Gravity (further – HMPG) is a very interesting choice as a modified gravity theory. In particular, HMPG completed some simple and classical tests in the Solar system [9], and also, it was shown that this type of gravity could describe an accelerated universe without dark energy (\(\Lambda \) term) [8]. HMPG and \(f({\mathscr {R}})\) are fourth-order theories of gravity and they are equivalent to two different scalar-tensor theories. Indeed, there exists a conformal transformation that connects the two theories.

2 Traversable wormholes in classical GR gravity

The static non-charged traversable wormhole proper line element, known as Morris–Thorne wormhole is given by the following expression [31]:

$$\begin{aligned} ds^2 = -e^{2\Omega (r)}dt^2+\frac{1}{1-\dfrac{b(r)}{r}}dr^2+r^2d\theta ^2 + r^2\sin ^2\theta d\phi ^2. \end{aligned}$$
(1)

Same as Eq. (1), but in Cartesian coordinates the line element becomes

$$\begin{aligned} ds^2 = -e^{2\Omega (x)}dt^2+\frac{1}{1-\dfrac{b(x)}{x}}dx^2+dy^2 + dz^2. \end{aligned}$$
(2)

Function \(\Omega (x)\) is the redshift function and b(x) is so-called wormhole shape function. For GR we have following Action Integral and gravitational Lagrangian \({\mathscr {L}}_g\)

$$\begin{aligned} {\mathscr {S}}_{\mathscr {R}}=\int _{\mathscr {M}}d^4x \sqrt{-g}\frac{1}{2}{\mathscr {L}}_g = \int _{\mathscr {M}}d^4x \sqrt{-g}\frac{1}{2}{\mathscr {R}}, \end{aligned}$$
(3)

where we have assumed for the gravitational constant \(\kappa =1\). In the latter Action Integral, g is metric tensor determinant: \(g=\det g_{\mu \nu }\), and \({\mathscr {R}}\) is the Ricci scalar of g. The Einstein field equations are (Einstein Field Equation or EFE):

$$\begin{aligned} G_{\mu \nu }+\Lambda g_{\mu \nu }=T_{\mu \nu }, \end{aligned}$$
(4)

where \(G_{\mu \nu }\) is the Einstein tensor, \(\Lambda \) is so-called lambda-term or dark energy (further – DE, in our case we consider universe without DE, so \(\Lambda =0\)), \(T_{\mu \nu }\) is stress-energy tensor for the additional matter source. The Einstein tensor is defined as follows

$$\begin{aligned} G_{\mu \nu } = R_{\mu \nu }-\frac{1}{2}{\mathscr {R}}g_{\mu \nu }. \end{aligned}$$
(5)

Here \(R_{\mu \nu }\) is the Ricci tensor, \({\mathscr {R}}\) is the Ricci scalar. So, now we can derive the non-zero components of the Einstein tensor for the line element (1):

$$\begin{aligned} G_{tt}= & {} -\frac{b'(r)}{r^2}, \end{aligned}$$
(6)
$$\begin{aligned} G_{rr}= & {} -\frac{\left( 1-\dfrac{b(r)}{r}\right) \left( 2 r b(r) \Omega '(r)+b(r)-2 r^2 \Omega '(r)\right) }{r^2 (r-b(r))}, \end{aligned}$$
(7)
$$\begin{aligned} G_{\theta \theta }= & {} \frac{\left( r \Omega '(r)+1\right) \left( -r b'(r)+2 r (r-b(r)) \Omega '(r)+b(r)\right) }{2 r^3} \nonumber \\&+\frac{(r-b(r)) \Omega ''(r)}{r}, \end{aligned}$$
(8)
$$\begin{aligned} G_{\phi \phi }= & {} \frac{\left( r \Omega '(r)+1\right) \left( -r b'(r)+2 r (r-b(r)) \Omega '(r)+b(r)\right) }{2 r^3}\nonumber \\&+\frac{(r-b(r)) \Omega ''(r)}{r}. \end{aligned}$$
(9)

In Eqs. (6), (7), (8) and (9) prime \(\Omega '(r)\), \(b'(r)\) means total derivative with respect to the independent variable r.

3 EFE’s for Morris–Thorne wormholes in different modified theories of gravity

3.1 \(f({\mathscr {R}}\)) modified gravity case

In \(f({\mathscr {R}})\) theory of gravity the EH Action Integral is modified as follows [6]:

$$\begin{aligned} \begin{aligned} {\mathscr {S}}_{\mathscr {R}}&=\frac{1}{2}\int _{{\mathscr {M}}} ({\mathscr {R}}+{\mathscr {L}}_M) \sqrt{-g} d^4x \\&\quad \Downarrow \\ {\mathscr {S}}_{f({\mathscr {R}})}&=\frac{1}{2}\int _{{\mathscr {M}}} [f({\mathscr {R}})+{\mathscr {L}}_M)]\sqrt{-g} d^4x, \end{aligned} \end{aligned}$$
(10)

where \({\mathscr {L}}_\text {M}\) is the matter Lagrangian. In the [47] we already modified general view of the EFE for symmetric metric tensor [48]:

$$\begin{aligned} \begin{aligned} G^{(0)}_{\mu \nu }\equiv R_{\mu \nu }-\frac{1}{2}g_{\mu \nu }{\mathscr {R}}&=\frac{T^\text {M}_{\mu \nu }}{f'({\mathscr {R}})}{+}g_{\mu \nu }\frac{[f({\mathscr {R}}){-}{\mathscr {R}}f'({\mathscr {R}})]}{2f'({\mathscr {R}})}\\&\quad + \frac{[\nabla _\mu \nabla _\nu f'({\mathscr {R}})-g_{\mu \nu }\Box f'({\mathscr {R}})]}{f'({\mathscr {R}})}. \end{aligned} \end{aligned}$$
(11)

Clearly \(f'({\mathscr {R}})\) is the derivative with respect to the Ricci scalar. Furthermore, the stress-energy tensor for the anisotropic fluid is [45]:

$$\begin{aligned} T^\text {M}_{\mu \nu }=\begin{pmatrix} -\rho &{}\quad 0 &{}\quad 0 &{}\quad 0\\ 0 &{}\quad p_r &{}\quad 0 &{}\quad 0\\ 0 &{}\quad 0 &{}\quad p_t &{}\quad 0\\ 0 &{}\quad 0 &{}\quad 0 &{}\quad p_t\\ \end{pmatrix}, \end{aligned}$$
(12)

where \(p_r\), \(p_t\) are the radial and tangential pressures respectively, and \(\rho \) is the energy-density. The parameter for the equation of state (further – EoS) defines the type of matter (\(\omega = \frac{p}{\rho }\)).

3.1.1 Violation of null energy condition (NEC)

In GR wormholes are supported by exotic matter, which involves a stress-energy tensor that violates the null energy condition (NEC) [27, 31, 51]. NEC violation in \(f({\mathscr {R}})\) is given by the following expression:

$$\begin{aligned} \begin{aligned} \rho ^{\text {eff}}+p_r^{\text {eff}}&=\frac{\rho +p_r}{f'({\mathscr {R}})}+\frac{1}{f'({\mathscr {R}})}\\&\bigg [(f'({\mathscr {R}}))''-(f'({\mathscr {R}}))'\frac{b'r-b}{2r^2(1-b/r)}\bigg ]\\&\quad \Downarrow \\ \rho ^{\text {eff}}+p_r^{\text {eff}}&=\frac{b'r-b}{r^3}<0. \end{aligned} \end{aligned}$$
(13)

Also, if a wormhole exists, it should obey the following inequalities [27]:

$$\begin{aligned} \frac{f'({\mathscr {R}})b'}{r^2} \ge 0 \Rightarrow \frac{b'}{r^2}<0 \Rightarrow f'({\mathscr {R}}) \le 0. \end{aligned}$$
(14)

Theorem of wormhole non-existence in \(f({\mathscr {R}})\) gravity is essentially the same the one derived before in [4]:

$$\begin{aligned} \frac{d f({\mathscr {R}})}{d{\mathscr {R}}} < 0. \end{aligned}$$
(15)

We continue by assuming specific functional form for the \(f({\mathscr {R}})\) function.

3.1.2 Derivation of the shape function

One of the necessary conditions for a wormhole to exist is that wormhole shape function must satisfy EFE’s in modified gravity. From Eq. (11) one could derive following EFE’s (in general \(f({\mathscr {R}})\) gravity with arbitrary choice of function) [27]:

$$\begin{aligned} \rho= & {} \frac{f_Rb'}{r^2}, \end{aligned}$$
(16)
$$\begin{aligned} p_r= & {} -\frac{b f_R}{r^3}+\frac{f'_R}{2r^2}(b'r-b)-f''_R\bigg (1-\frac{b}{r}\bigg ), \end{aligned}$$
(17)
$$\begin{aligned} p_t= & {} -\frac{f'_R}{r}\bigg (1-\frac{b}{r}\bigg )+\frac{f_R}{2r^3}(b-b'r). \end{aligned}$$
(18)

Here we already considered wormhole solution without tidal forces, and thus with constant redshift factor. We have chosen this case, because the tidal gravitational forces experienced by a traveler must be bearably small (negligible) [31]. And, therefore with applying equation of state \(p_t=\omega \rho \) we have following equation [27]:

$$\begin{aligned} f'_R\bigg (1-\frac{b}{r}\bigg )-\frac{f_R}{2r^2}[b-b'r(1+2\omega )] = 0. \end{aligned}$$
(19)

Hence, from equation above we could find suitable shape function which satisfies EFE’s for any kind of \(f({\mathscr {R}})\) gravity.

3.1.3 Quadratic gravity \(f({\mathscr {R}})={\mathscr {R}}+\alpha {\mathscr {R}}^2\)

Firstly, we could rewrite Eq. (13) in terms of quadratic MOG:

$$\begin{aligned} \begin{aligned}&\left( 1-\frac{b(r)}{r}\right) \left( \frac{4 \alpha b''(r)}{r^2}-\frac{8 \alpha b'(r)}{r^3}\right) \\&\quad -\frac{\left( \frac{4 \alpha b'(r)}{r^2}+1\right) \left( b(r)-r (2 \omega +1) b'(r)\right) }{2 r^2}=0. \end{aligned} \end{aligned}$$
(20)

We solved this equation numerically (there is no possibility to solve this equation algebraically) with initial conditions \(b(10^{-2})=10^{-3}\), \(b'(10^{-2})=2\times 10^{-4}\) [22].

As well, for our MOG shape function, we have following flaring-out condition [23]:

$$\begin{aligned} \begin{aligned} \frac{(b-b'r)}{b^2}>0 \Rightarrow b'(r_0)<1. \end{aligned} \end{aligned}$$
(21)
Fig. 2
figure 2

Numerical solution for (20) and Flaring-out condition validation/violation for positive \(\alpha \) parameter in quadratic gravity

We numerically solved Eqs. (20) and (21) on the Figure. 2. The flaring-out condition was solved for only positive values of \(\alpha \), because, as it turned out, with \(\alpha <0\) we have that flaring-out condition is violated generally. Thus, in quadratic MOG for positive values of MOG parameter, we have physically acceptable shape function, that satisfies EFE’s. Finally, for almost all positive alpha values, the flaring-out condition was satisfied with EoS parameter \(\omega = 1\) (stiff fluid, presented by [54]) Therefore, we could proceed to the NEC, WEC, and SEC conditions derivation for this MOG, shape function.

3.1.4 Quadratic gravity energy conditions

Because of the previously stated reasons (validation of flare-out condition for every \(\alpha \ge 0\) at the throat), we chose the case with the stiff fluid (as numerical analysis showed, there is no significant differences of energy conditions with different values of EoS parameter in the limit \(0<\omega \le 1\)). Firstly, we could present energy conditions, that we consider in this paper [1, 44]:

  • Null Energy Condition (NEC): \(\rho +p_r \ge 0\wedge \rho +p_t \ge 0\)

  • Weak Energy Condition (WEC): \(\rho -p_r \ge 0 \wedge \rho +p_t \ge 0\)

  • Strong Energy Condition (SEC): \(\rho + p_r + 2p_t \ge 0\)

NEC is a minimal requirement of WEC and SEC conditions and must be obeyed always (if NEC is violated, so-called exotic matter or in some cases phantom fluid will appear [41]).

Fig. 3
figure 3

NEC, WEC and SEC conditions for quadratic gravity with \(\omega =1\)

We showed numerical solutions for NEC, WEC and SEC energy conditions at Figure. 3. As one may notice, generally NEC is validated for each pressure if \(\alpha <1\), as well WEC is validated for both pressure types for any \(\alpha >0\). Finally, SEC is violated.

3.1.5 Power law gravity \(f({\mathscr {R}})=f_0{\mathscr {R}}^n\) for \(n>1\)

This is another \(f({\mathscr {R}})\) example of MOG theory, which was described by [10, 34]. For this type of gravity we have following EH action:

$$\begin{aligned} {\mathscr {S}}_{f({\mathscr {R}})}=\frac{1}{2}\int _{\mathscr {M}} [f_0{\mathscr {R}}^n]\sqrt{-g}d^4x, \end{aligned}$$
(22)

where \(f_0\) is a constant to give correct dimensions to the action and n is the slope parameter [28]. For MOG of this kind we have following form of Eq. (19):

$$\begin{aligned} \begin{aligned} f_0 2^{n-1} (n-1) n \left( 1-\frac{b(r)}{r}\right) \left( \frac{b'(r)}{r^2}\right) ^{n-2} \left( \frac{b''(r)}{r^2}-\frac{2 b'(r)}{r^3}\right) \\ -\frac{f_0 2^{n-2} n \left( \frac{b'(r)}{r^2}\right) ^{n-1} \left( b(r)-r (2 \omega +1) b'(r)\right) }{r^2}=0. \end{aligned} \end{aligned}$$
(23)
Fig. 4
figure 4

Numerical solution for Eq. (23) and Flaring-out condition validation/violation for varying n and \(\omega \) in power-law gravity

On the Fig. 4 we once again plotted numerical solutions for Eqs. (23) and (13) with varying n/constant \(\omega \) and varying \(\omega \)/constant n. As it turned out, the shape function does not depends on the \(f_0\) parameter. Now that we have decided on the type of shape function for our power-law MOG, we can begin to study the energy conditions of the Morris–Thorne traversable wormhole in the MOG of our consideration.

Fig. 5
figure 5

NEC, WEC and SEC conditions for power-law gravity with \(\omega =1\) and \(n=3\)

In turn, in the Fig. 5 we illustrated null, weak, and strong energy conditions for power-law gravity. We plotted only one case with \(\omega =1\) and \(n=3\), because we found that if we will vary these parameters, nothing changes much. Also, it is interesting that our power-law numerical solutions for energy conditions are very similar to those, that we obtained for quadratic MOG (see Fig. 3). In relation, NEC generally is also validated for both pressures, but now with \(\alpha <2\), for \(\alpha <4\) WEC is validated for both pressure types and finally SEC is violated.

3.1.6 Logarithmic corrected \(f({\mathscr {R}})\) gravity

Logarithmic corrected \(f({\mathscr {R}})\) gravity was introduced in [40]:

$$\begin{aligned} f({\mathscr {R}})={\mathscr {R}}+\alpha {\mathscr {R}}^2+\beta {\mathscr {R}}^2\log \beta {\mathscr {R}}, \end{aligned}$$
(24)

where \(\beta >0\) and \(\alpha >0\). This type of gravity can describe expanding universe without dark energy [34]. The modified EH Action Integral is defined in this way:

$$\begin{aligned} S_{f({\mathscr {R}})}=\int _{\mathscr {M}}d^4x\sqrt{-g}\frac{1}{2}\bigg ({\mathscr {R}}+\alpha {\mathscr {R}}^2+\beta {\mathscr {R}}^2\log \beta {\mathscr {R}}\bigg ), \end{aligned}$$
(25)

Now, with given MOG form we as usual could rewrite Eq. (19) as follows:

$$\begin{aligned}&2 (2 \omega +1) b'(r)^2 \bigg (2 \alpha +2 \beta \log \bigg (\frac{2 \beta b'(r)}{r^2}\bigg )+\beta \bigg )\nonumber \\&\quad +b'(r) \bigg (-8 (2 \alpha +3 \beta )-16 \beta \log \bigg (\frac{2 \beta b'(r)}{r^2}\bigg )+r^2 (2 \omega +1)\bigg )\nonumber \\&\quad +4 r b''(r) \bigg (2 \alpha +2 \beta \log \bigg (\frac{2 \beta b'(r)}{r^2}\bigg )+3 \beta \bigg )\nonumber \\&\quad +b(r) \Bigg (\frac{2 b'(r) \bigg (6 \alpha +6 \beta \log \bigg (\frac{2 \beta b'(r)}{r^2}\bigg )+11 \beta \bigg )}{r}\nonumber \\&\quad -4 b''(r) \bigg (2 \alpha +2 \beta \log \bigg (\frac{2 \beta b'(r)}{r^2}\bigg )+3 \beta \bigg )-r\Bigg )=0. \end{aligned}$$
(26)
Fig. 6
figure 6

Numerical solution for Eq. (26) and Flaring-out condition validation/violation for varying \(\alpha \) and \(\omega \) (we assumed constant \(\beta \) for simplicity) in log-corrected gravity

Consequently, on the Fig. 6 we have numerical representation of Eq. (26) and of flare-out condition at the throat. As numerical analysis showed, \(b'(r_0)<1\) is obeyed for \(0<\omega \le 1\wedge \alpha >4\) if we assume constant \(\beta =1\).

Fig. 7
figure 7

NEC, WEC and SEC conditions for log-corrected gravity with \(\omega =\beta =1\)

Routinely, for this shape function we placed NEC, WEC and SEC numerical solutions at the Fig. 7. As we found, NEC is violated for radial pressure and validated for tangential, WEC is vice versa validated for radial pressure case and violated for tangential one. Just as with the previous MOG theories, SEC is generally violated.

3.2 Hybrid metric-Palatini gravity

The second family of modified theories of our consideration is the HMPG theory, for which the Action Integral is given by [8, 18, 37]:

$$\begin{aligned} {\mathscr {S}}_{f(\mathscr {\hat{R}})}=\frac{1}{2}\int _{\mathscr {M}}d^4x \sqrt{-g}[{\mathscr {R}}+f(\mathscr {\hat{R}})], \end{aligned}$$
(27)

where \(\mathscr {\hat{R}}\) is the Palatini scalar, constructed from Palatini curvature tensor, which reads [25]:

$$\begin{aligned} \hat{R}_{\mu \nu } = \hat{R}^\alpha _{\mu \alpha \nu }= {\hat{\Gamma }}^\alpha _{\mu \nu ,\alpha }-{\hat{\Gamma }}^\alpha _{\mu \alpha ,\mu } + {\hat{\Gamma }}^\alpha _{\alpha \lambda }{\hat{\Gamma }}^\lambda _{\mu \nu }-{\hat{\Gamma }}^\alpha _{\mu \lambda }{\hat{\Gamma }}^\lambda _{\alpha \nu }. \end{aligned}$$
(28)

\({\hat{\Gamma }}^\alpha _{\mu \nu }\) is the Levi-Civita connection for metric which is conformal to our wormhole background metric (\(h_{\mu \nu }=\phi g_{\mu \nu }\))[3]. We could also rewrite Eq. (28) in form [12]:

$$\begin{aligned} \hat{R}_{\mu \nu } = R_{\mu \nu }+\frac{3}{2\phi ^2}\partial _\mu \phi \partial _\nu \phi -\frac{1}{\phi }\nabla _\mu \nabla _\nu \phi -\frac{1}{2\phi }g_{\mu \nu }\Box \phi . \end{aligned}$$
(29)

Thus for the Palatini scalar we find that

$$\begin{aligned} \mathscr {\hat{R}} = {\mathscr {R}}+\frac{3}{2\phi ^2}\partial _\mu \phi \partial ^\mu \phi - \frac{3}{\phi }\Box \phi . \end{aligned}$$
(30)

But, we have as well scalar-tensor representation of hybrid metric-Palatini gravity [19]:

$$\begin{aligned} {\mathscr {S}}_{f(\mathscr {\hat{R}})}=\frac{1}{2}\int _{\mathscr {M}}d^4x\sqrt{-g}\bigg [(1+\phi ){\mathscr {R}}+\frac{3}{2\phi }\partial _\mu \phi \partial ^\mu \phi -V(\phi )\bigg ], \end{aligned}$$
(31)

where \(\phi \) is scalar field and \(V(\phi )\) is scalar potential (in metric-Palatini gravity scalar field in dynamic [35]). By varying the action, we could obtain following EFE form [19, 25]:

$$\begin{aligned} G_{\mu \nu }= & {} \frac{1}{1+\phi }T_{\mu \nu }+T^{(\phi )}_{\mu \nu }, \end{aligned}$$
(32)
$$\begin{aligned} T^{(\phi )}_{\mu \nu }= & {} \frac{1}{1+\phi }\Bigg [\nabla _\mu \nabla _\nu \phi -\frac{3}{2\phi }\nabla _\mu \phi \nabla _\nu \phi \nonumber \\&+ \Bigg (\frac{3}{4\phi }\nabla _\lambda \phi \nabla _\lambda \phi - \Box \phi -\frac{1}{2}V(\phi )\Bigg )\Bigg ], \end{aligned}$$
(33)

where \(V(\phi )\) is [3]:

$$\begin{aligned} V(\phi ) = \mathscr {\hat{R}}\phi -f(\mathscr {\hat{R}}), \end{aligned}$$
(34)

and scalar field is [12]:

$$\begin{aligned} \phi = \tan ^2\bigg (\sqrt{\frac{3}{8}\bar{\phi }}\bigg ). \end{aligned}$$
(35)

If \(\phi \rightarrow \infty \), then:

$$\begin{aligned} \bar{\phi } = \sqrt{\frac{8}{3}}\bigg [(-1)^k\frac{\pi }{2}+2k\pi \bigg ], \qquad k=0,1,2,3\ldots \end{aligned}$$
(36)

On the other hand, when \(\phi \rightarrow 0\),

$$\begin{aligned} \bar{\phi } = \sqrt{\frac{8}{3}}k\pi . \end{aligned}$$
(37)

If (from Eqs. (35) and (36)) \(\phi \) is independent of \(x^\mu \), then \(\mathscr {\hat{R}}={\mathscr {R}}\). As we did for \(f({\mathscr {R}})\) gravity, we could derive energy density and pressures from EFE’s:

$$\begin{aligned} \rho= & {} \frac{-2 b'(r)-f(\hat{{\mathscr {R}}})r^2}{2 r^2}, \end{aligned}$$
(38)
$$\begin{aligned} p_r= & {} \bigg [2 (b(r)-r) \cot ^2(2 \pi k) b'(r)-b(r) (f(\hat{{\mathscr {R}}})r^2 \nonumber \\&+2 \cot ^2(2 \pi k)+2)+f(\hat{{\mathscr {R}}})r^3\bigg ]\bigg /\bigg [2 r^2 (r-b(r))\bigg ], \end{aligned}$$
(39)
$$\begin{aligned} p_t= & {} \frac{1}{2} \Bigg (\bigg [\csc ^2(2 \pi k) \bigg (r b(r)-b'(r) (\cos (4 \pi k) \nonumber \\&+r^2+1\bigg )\bigg )\bigg ]\bigg /[r^2]+f(\hat{{\mathscr {R}}})\Bigg ). \end{aligned}$$
(40)

Then, EoS \(p_t/\rho =\omega \) takes form:

$$\begin{aligned}&\bigg [\csc ^2(2 \pi k) \left( 2 b'(r) \left( \cos (4 \pi k)+r^2+1\right) -2 r b(r)\right. \nonumber \\&\quad \left. +f(\hat{{\mathscr {R}}}) r^2 (\cos (4 \pi k)-1)\right) \bigg ]\bigg /\bigg [4 b'(r)+2 f(\hat{{\mathscr {R}}}) r^2\bigg ] = \omega . \nonumber \\ \end{aligned}$$
(41)

By solving this equation, we could obtain shape function b(r) in the physically acceptable form (i.e. shape function that satisfies EFE’s)

3.2.1 Exponential \(f(\mathscr {\hat{R}})=\zeta \bigg (1+e^{-\frac{\mathscr {\hat{R}}}{\varPhi }}\bigg )\)

In this study we consider the following Palatini-scalar function [17]:

$$\begin{aligned} f(\mathscr {\hat{R}}) = \zeta \bigg (1+e^{-\frac{\mathscr {\hat{R}}}{\varPhi }}\bigg ). \end{aligned}$$
(42)

Therefore, we could rewrite EoS (41), that describe shape function, which satisfies EFE’s as follows:

$$\begin{aligned}&\bigg [\csc ^2(2 \pi k) \left( r b(r)-b'(r) \left( \cos (4 \pi k)+r^2+1\right) \right) \nonumber \\&\quad +\zeta r^2 \left( e^{-\frac{2 b'(r)}{r^2 \varPhi }}+1\right) \bigg ]\bigg /\bigg [\zeta r^2 \left( e^{-\frac{2 b'(r)}{r^2 \varPhi }}+1\right) +2 b'(r)\bigg ] = \omega .\nonumber \\ \end{aligned}$$
(43)
Fig. 8
figure 8

Shape function and flare-out condition at the throat for HMPG gravity with \(\omega =-1\) and \(k=\pi /2\)

Hence, on the Fig. 8 we located the numerical solution for the equation above and proof that in this MOG for our shape function flare-our condition is satisfied. As we saw from the numerical analysis, shape function does not depends on the \(\zeta \) MOG free parameter. Also, it is necessary to note that we have only one physically acceptable solution with EoS parameter \(\omega =-1\) (dark energy like fluid).

Fig. 9
figure 9

NEC, WEC and SEC conditions for hybrid metric-Palatini gravity with \(\omega =-1\), \(\varPhi =1\) and \(k=\pi /2\)

Finally, in the Fig. 9 we illustrated energy conditions for HMPG gravity. As one may notice, NEC is violated for both pressure types, and thus there is always present exotic matter at the throat. WEC condition could be satisfied by assuming that \(\zeta <1\), SEC is validated for every \(\zeta >0\).

4 Quantization of exotic matter, that violate NEC condition

Volume integral quantifier (further – just VIQ) could help us with the derivation of exact exotic matter volume. With VIQ we have the opportunity to understand with which values of the MOG parameters the volume of matter violating the null energy condition is the smallest in the case where exotic matter is present [32].

Volume integral quantifier is given by [42]:

$$\begin{aligned} \varPhi= & {} \int ^\infty _{r_0}\int ^\pi _0\int ^{2\pi }_0[\rho +p_r]\sqrt{-g}drd\theta d\phi \nonumber \\&\Downarrow \nonumber \\&\oint [\rho +p_r]dV = 2\int ^\infty _{r_0}[\rho +p_r]4\pi r^2 dr. \end{aligned}$$
(44)
Fig. 10
figure 10

Volume integral quantifier for quadratic, power-law, log-corrected and HMPG gravities. For the \(f({\mathscr {R}})\) MOG’s we considered case with \(\omega =1\) and for HMPG one with \(\omega =-1\). Also, for both kinds of theories we assumed \(r_0=1\)

From Fig. 10 obviously for any choice of MOG parameter if \(r\rightarrow \infty \) then \(\varPhi \rightarrow 0\). Also, for any \(\alpha>0\wedge f_0>0\wedge \zeta >0\) we have exotic matter at the throat, but we could minimize its amount in the first two MOG’s of our consideration if we set \(\alpha \rightarrow 0\) and \(f_0\rightarrow 0\). For log-corrected gravity \(\varPhi \) is minimized at \(\alpha =5\), for hybrid metric-Palatini gravity we couldn’t minimize exotic matter contribution because \(\varPhi \) is independent of \(\zeta \) MOG parameter.

5 Wormhole stability in modified theories of gravity

Wormhole stability conditions can be examined by employing an equilibrium condition obtained from the Tolman–Oppenheimer–Volkov equation for non-tidal traversable wormhole: [36, 38, 39, 50]:

$$\begin{aligned} \begin{aligned} \frac{dp_r}{dr}+\underbrace{\Omega (r)'(\rho +p_r)}_{\mathrm{if}}\ \Omega&=0,0+\frac{2}{r}(p_r-p_t) + F_{\mathrm {ex}}=0\\&\quad \Downarrow \\&\frac{dp_r}{dr}+\frac{2}{r}(p_r-p_t) + F_{\mathrm {ex}} =0. \end{aligned} \end{aligned}$$
(45)

One may see that in our modified TOV (further – MTOV) present extra force \( F_{\mathrm {ex}}\), which exists to hold WH stable [13, 46], even if stress-energy tensor is not conserved (in considered theories of gravity generally \(\nabla ^\mu T_{\mu \nu }\ne 0\) [29, 30, 43]). If the wormhole is stable, the MTOV conditions must be satisfied. Hence, for each model of our analysis it follows.

Fig. 11
figure 11

MTOV forces example in quadratic gravity with \(\alpha =\omega =1\)

Routinely, on the Fig. 11 we have located the example of forces, that present in MTOV (45) for quadratic gravity case with \(\alpha =\omega =1\).

5.1 \(f({\mathscr {R}})={\mathscr {R}}+\alpha {\mathscr {R}}^2\) gravity

To satisfy MTOV, external force must look like:

$$\begin{aligned} \begin{aligned} F_{\mathrm {ex}}&=\bigg [2 \alpha \bigg (2 r^3 (r-b(r)) b^{(4)}(r)-18 r b'(r)^2 +r (r (11 b(r)\\&\quad -8 r) b^{(3)}(r)-b''(r) \left( r \left( r b''(r)-16\right) +26 b(r)\right) ) +b'(r)\\&\quad \times \left( r \left( -3 r^2 b^{(3)}(r)+14 r b''(r)-16\right) +30 b(r)\right) \bigg )\bigg ]\bigg /[r^6]. \end{aligned} \end{aligned}$$
(46)
Fig. 12
figure 12

MTOV extra force with varying \(\alpha \) for \(\omega =1\)

On the Fig. 12 we numerically solved Eq. (46). Judging by data from the figure above, we could say that as \(\alpha \rightarrow 0\Rightarrow F_{\mathrm {ex}}\rightarrow 0\) (GR restored). On other hand, for non-zero MOG parameter to keep wormhole stable there is always must be present extra force.

5.2 \(f({\mathscr {R}})=f_0{\mathscr {R}}^n\) gravity

MTOV extra force for power law gravity have following form:

$$\begin{aligned}&F_{\mathrm {ex}} = \Bigg [f_0 2^{n-2} (n-1) n \left( \frac{b'(r)}{r^2}\right) ^{n+1} (2 (n-3) (n-2) \nonumber \\&\quad \times r^3 (r-b(r)) b''(r)^3+6 (1-2 n) r b'(r)^4 \nonumber \\&\quad -b'(r)^3 \left( r \left( r \left( 3 r b^{(3)}(r)+2 (5-6 n) b''(r)\right) +8 n (2 n-3)\right) \right. \nonumber \\&\quad \left. +2 \left( -8 n^2+6 n+5\right) b(r)\right) +(n-2) r^2 b'(r) b''(r) \nonumber \\&\quad \times \left( 6 r (r-b(r)) b^{(3)}(r)+((12 n-13) b(r)+4 (4-3 n) r) b''(r)\right) \nonumber \\&\quad +r b'(r)^2 \left( (5-3 n) r^2 b''(r)^2+2 ((4 (5-3 n) n-5) b(r)\right. \nonumber \\&\quad \left. +2 (n (6 n-13)+6) r) b''(r)+r (2 r (r-b(r)) b^{(4)}(r)\right. \nonumber \\&\quad \left. +b^{(3)}(r) ((12 n-13) b(r)+4 (4-3 n) r)\right) ))\Bigg ]/[b'(r)^5]. \end{aligned}$$
(47)
Fig. 13
figure 13

MTOV extra force with varying \(f_0\) for \(\omega =1\)

As usual, we depicted the solution for Eq. (47) on the Fig. 13.

5.3 \(f({\mathscr {R}})={\mathscr {R}}+\alpha {\mathscr {R}}^2+\beta {\mathscr {R}}^2\ln \beta {\mathscr {R}}\) gravity

Analytical solution for MTOV extra force in log-corrected gravity:

$$\begin{aligned}&F_{\mathrm {ex}} = \bigg [4 \beta r^3 (b(r)-r) b''(r)^3-6 r b'(r)^4 \nonumber \\&\quad \times \bigg (6 \alpha +6 \beta \log \left( \frac{2 \beta b'(r)}{r^2}\right) +13 \beta \bigg )+b'(r)^3 (2 b(r) \nonumber \\&\quad \times \bigg (30 \alpha +30 \beta \log \left( \frac{2 \beta b'(r)}{r^2}\right) +97 \beta \bigg ) \nonumber \\&\quad +r \bigg (r \bigg (b''(r) \bigg (28 \alpha +28 \beta \log \bigg (\frac{2 \beta b'(r)}{r^2}\bigg )+66 \beta \bigg ) \nonumber \\&\quad -3 r b^{(3)}(r)\bigg (2 \alpha +2 \beta \log \bigg (\frac{2 \beta b'(r)}{r^2}\bigg )+3 \beta \bigg )\bigg )\nonumber \\&\quad -32 \bigg (\alpha +\beta \log \bigg (\frac{2 \beta b'(r)}{r^2}\bigg )+4 \beta \bigg )\bigg )\bigg )+2 \beta r^2 b'(r) b''(r)\nonumber \\&\quad \times (6 r (r-b(r)) b^{(3)}(r)+(11 b(r)-8 r) b''(r))+r b'(r)^2\nonumber \\&\quad \times \bigg (-r^2 b''(r)^2 (2 \alpha +2 \beta \log \bigg (\frac{2 \beta b'(r)}{r^2}\bigg )+9 \beta \bigg )+2 b''(r)\nonumber \\&\quad \times \bigg (2 \beta (8 r-13 b(r)) \log \bigg (\frac{2 \beta b'(r)}{r^2}\bigg )-26 \alpha b(r)\nonumber \\&\quad -95 \beta b(r)+16 \alpha r+68 \beta r\bigg )+r \bigg (2 r (r-b(r)) b^{(4)}(r) \nonumber \\&\quad \times \bigg (2 \alpha +2 \beta \log \bigg (\frac{2 \beta b'(r)}{r^2}\bigg )+3 \beta \bigg )\nonumber \\&\quad +b^{(3)}(r) \bigg (2 \beta (11 b(r) -8 r\bigg ) \log \bigg (\frac{2 \beta b'(r)}{r^2}\bigg )\nonumber \\&\quad +(22 \alpha +57 \beta ) b(r)-16 r (\alpha +3 \beta )\bigg )\bigg )\bigg )\bigg ]\bigg /\bigg [r^6 b'(r)^2\bigg ].\nonumber \\ \end{aligned}$$
(48)
Fig. 14
figure 14

MTOV extra force with varying \(\alpha \) for \(\beta =\omega =1\)

Therefore, as we did for other MOG’s, we placed solution of Eq. (48) on the Fig. 14. Case with zero extra force (GR) could be obtained by setting \(\alpha =\beta =0\).

5.4 Metric-Palatini exponential gravity

For this MOG we have following \(F_{\mathrm {ex}}\):

$$\begin{aligned}&F_{\mathrm {ex}} = \Bigg [\bigg [\csc ^2(2 \pi k) (2 r b(r) (-2 r \cos ^2(2 \pi k) b''(r)+2 b'(r)\nonumber \\&\quad \times (\cos (4 \pi k)+r^2+1)+r^2-1)+b(r)^2 (r ((\cos (4 \pi k)+1) b''(r)\nonumber \\&\quad -4 r)-2 b'(r) (\cos (4 \pi k)+r^2+1))+r^2 (r (\cos (4 \pi k)+1) \nonumber \\&\quad \times b''(r)-2 b'(r) (\cos (4 \pi k)+r^2))+2 r b(r)^3)\bigg ]\bigg /\bigg [(r-b(r))^2\bigg ]\nonumber \\&\quad +\frac{2 \zeta e^{-\frac{2 b'(r)}{r^2 \varPhi }} (r b''(r)-2 b'(r))}{\varPhi }\Bigg ]\Bigg /\Bigg [2 r^3\Bigg ]. \end{aligned}$$
(49)
Fig. 15
figure 15

MTOV extra force with varying \(\zeta \) for \(\varPhi =1\) and \(\omega =-1\)

Finally, we numerically solved Eq. (49) for last MOG of our consideration (HMPG) at the Fig. 15.

6 Conclusions

We presented Morris–Thorne traversable wormhole solutions for different modified gravity theories, such as: \(f({\mathscr {R}}) = {\mathscr {R}}+a{\mathscr {R}}^2\), \(f({\mathscr {R}}) = f_0{\mathscr {R}}^n\), \(f({\mathscr {R}})={\mathscr {R}}+\alpha {\mathscr {R}}^2+\beta {\mathscr {R}}^2\ln \beta {\mathscr {R}}\) and hybrid metric-Palatini gravity \(f(\mathscr {\hat{R}})\). For each kind of modified gravity we derived suitable shape function that satisfies Einstein Field Equation by applying EoS \(p_t=\omega \rho \). Numerical solutions for b(r) and \(b'(r_0)\) are represented at Figs. 2, 4, 6 and 8.

We probed these models via numerical solutions of the null energy, weak energy and strong energy conditions, for \(f({\mathscr {R}})\) family of gravity theories results are presented at Figs. 3, 5 and 7, for metric-Palatini gravity at Fig. 9.

Moreover, we found a volume integral quantifier. The latter was used to construct plots with a volume of matter that violates NEC condition, i.e. exotic matter, and to obtain some parameter values and conditions, which could reduce the amount of exotic matter near the wormhole throat. The results of matter quantifying can be recognized in Fig. 10. As well, we probed the stability of the non-tidal wormholes in the modified gravities by the modified Tolman–Oppenheimer–Volkov equation (equilibrium). Furthermore, we found the contribution of the extra force that arises because of the non-continuity of stress-energy tensor for each MOG w.r.t. free parameters. For the graphical representation of MTOV extra force solutions, see Figs. 12, 13, 14 and 15. More information about wormhole stability and suitable parameter values could be found in Sect. 5.

This study contributes to the subject of the existence of wormhole solutions in modified theories of gravity. For the two fourth-order theories of our consideration, we found that the HMPG provides wormhole solutions with a fluid source with value for the EoS close to the cosmological constant. In a future study, we plan to investigate the relation of these solutions under the action of the conformal transformation which relates the two theories.