Implementing and Experimenting with the Calabi-Wilf algorithm

for random selection of a subspace over a finite field

Shalosh B. EKHAD and Doron ZEILBERGER

In fond memory of Eugenio Calabi (May 11, 1923 - September 25, 2023), and
Herbert Saul Wilf (June 13, 1931 - January 7, 2012)

Abstract: In their beautiful note, “On the Sequential and Random Selection of
Subspaces over a Finite Field”, geometrical giant Eugenio Calabi and
combinatorial giant Herbert Wilf proposed an elegant algorithm to do what is
promised in their title. In the present note, written forty six years later,
we describe a Maple package, written by the second author, that implements
their beautiful algorithm, and then go on to report numerous experiments,
performed by the first author, that demonstrate the efficiency and reliability
of the Calabi-Wilf algorithm.

Memory Lane (by DZ)

One of my numerous temporary positions was during the academic year 1982-1983.
Herbert Wilf, who later became my close collaborator, got me a one-year
visiting position at the University of Pennsylvania. I occupied the office of
David Harbater, who was on sabbatical that year, and it so happened that it
was close to both Wilf’s and Calabi’s offices. I was in awe of both of them,
since they were already legendary way back then. One day Herb mentioned to me
that, while Calabi’s expertise was almost diametrically opposite to his, they
did have a non-empty intersection, and coauthored a cute little note [CW] that
applies the Wilf methodology [W] for random selection of combinatorial
objects, to something that besides being combinatorial (after all ‘finite
fields’ are finite!) was also a geometrical object, that Calabi can relate to:
namely subspaces of $GF(q)^{n}$.

Reminder about the Wilf approach for Random Generation of Combinatorial
objects

The binomial coefficient,

${{n}\choose{k}}={n!\over k!(n-k)!}\quad,$

famously counts the number of subsets with $k$ elements of a ‘universal’ set
of $n$ elements,$\\{1,\dots,n\\}$. It famously satisfies the Pascal-Chu
recurrence

${{n}\choose{k}}={{n-1}\choose{k-1}}\,+\,{{n-1}\choose{k}}\quad,$

that breaks up all $k$-element subsets of $\\{1,\dots,n\\}$ into those that do
contain $n$, whose number is the first term on the right, and those that do
not contain $n$, whose number is the second term on the right.

It follows that the fraction of $k$-element subsets of $\\{1,\dots,n\\}$ that
do contain $n$ is

${{{n-1}\choose{k-1}}\over{{n}\choose{k}}}={(n-1)!/((k-1)!(n-k)!)\over(n!/k!(n-k)!)}={k\over
n}\quad,$

hence the probability that a random $k$-subset of $\\{1,\dots,n\\}$ does
contain $n$ is ${k\over n}$, while the probability that it does not is
$1-{k\over n}$.

This simple observation lead Wilf [W] to suggest the following recursive
algorithm for the random generation, uniformly at random, of a $k$-subset of
$\\{1,\dots,n\\}$.

RandomSubset(n,k): if $k<0$ or $k>n$ return FAIL. If $k=0$ then RETURN the
empty set. Otherwise toss a loaded coin whose probability of Heads is $k/n$.
If it lends on Heads, output RandomSubset(n-1,k-1)$\cup\\{n\\}$, else output
RandomSubset(n-1,k) .

This philosophy was used on many other combinatorial objects that posses
recursive structures, and nicely described, and implemented (with FORTRAN
source-code included!), in the Nijenhuis-Wilf classic [NW]. However, the
Calabi-Wilf algorithm did not make it to that book. It was exposited in more
detail by Igor Pak, in his insightful ‘big picture’ paper [P], written in his
inimitable engaging style.

The Calabi-Wilf algorithm

The q-binomial coefficient,

$\left[\matrix{n\cr k}\right]={[n]!\over[k]![n-k]!}\quad,$

(where $[n]!:=(1-q)\cdots(1-q^{n})/(1-q)^{n}$ is the $q$-analog of $n!$), not-
quite-as-famously, counts the number of $k$-dimensional subspaces of
$GF(q)^{n}$. Each such subspace has many possible bases, of course, but only
one in row-echelon form, so in order to count $k$-dimensional subspaces one
should count $k\times n$ matrices over $GF(q)$ in row-echelon form.

The right side of the q-Pascal-Chu recurrence

$\left[\matrix{n\cr k}\right]\,=\,\left[\matrix{n-1\cr
k-1}\right]\,+\,q^{k}\left[\matrix{n-1\cr k}\right]\quad,$

may be thought of (per [CW]) as counting complementary subsets of the totality
of $k\times n$ matrices in echelon form:

$\bullet$ The first term counts those $k\times n$ matrices, $B=(b_{ij})_{1\leq
i\leq k,1\leq j\leq n}$ for which $b_{11}=1$, $b_{i1}=0$ ($1<i\leq k$),
$b_{1j}=0$ ($1<j\leq n)$, and for which necessarily the remaining
$(k-1)\times(n-1)$ array is a $(k-1)\times(n-1)$ matrix in row-echelon form.

$\bullet$ The second term counts the other basis matrices $B$, i.e. those
whose first column is an arbitrary vector of length $k$, and for which the
remaining $k\times(n-1)$ array is a basis matrix for a $k$ dimensional
subspace of $GF(q)^{n-1}$.

It follows that the fraction of matrices of the first kind is

$\left[\matrix{n-1\cr k-1}\right]/\left[\matrix{n\cr k}\right]={q^{k}-1\over
q^{n}-1}\quad,$

and it follows, in turn, that the probability that a random $k\times n$ matrix
in row-echelon form is of the first kind is ${q^{k}-1\over q^{n}-1}$.

This observation lead Calabi and Wilf [CW] to suggest the following recursive
algorithm for the random generation, uniformly at random, of a $k\times n$
matrix over $GF(q)$ in row-echelon form.

RandomSubspace(n,k): if $k<0$ or $k>n$ return FAIL. If $k=0$ then RETURN the
empty matrix (corresponding to the null subspace). Otherwise toss a loaded
coin whose probability of Heads is $(q^{k}-1)/(q^{n}-1)$. If it lends on
Heads, generate RandomSubspace(n-1,k-1), and convert it to a $k\times n$
matrix by putting on top the row vector $10^{n-1}$, and prefixing the column
vector $(10^{k-1})^{T}$ (where of course the sole $1$ is shared by both).

Else generate RandomSubspace(n-1,k). Also generate a random column vector in
$GF(q)$, of length $k$, and prefix it at the front.

Maple Implementation

The Calabi-Wilf algorithm (and much more!) is implemented in the Maple
package, CalabiWilf.txt, accompanying the present paper:

https://sites.math.rutgers.edu/~zeilberg/tokhniot/CalabiWilf.txt .

It is procedure Rqnk(q,n,k). For example, typing

Rqnk(7,10,5);

would generate, uniformly-at-random, one of the

$((7^{10}-1)(7^{9}-1)\cdots(7^{6}-1))/((7^{5}-1)(7^{4}-1)\cdots(7-1))=1602592475815614015216$

$5\times 10$ row-echelon matrices over $GF(7)$. Here is one of them:

$\left(\matrix{5&6&3&2&5&1&0&0&0&0\cr 6&1&0&0&6&0&1&0&0&0\cr
5&3&4&2&0&0&0&1&0&0\cr 0&5&1&6&6&0&0&0&1&0\cr
5&6&2&1&5&0&0&0&0&1}\right)\quad.$

Experimentation

$\bullet$ If you want to see an article about estimating the average,
variance, skewness, and kurtosis, for the number of occurrences of $1$s, in
row-echelon matrices of dimension $k\times 2k$ over GF(2), for $k$ from $50$
to $100$, by simulating, $1000$ times, look at the output file

http://www.math.rutgers.edu/~zeilberg/tokhniot/oCalabiWilf1.txt .

$\bullet$ If you want to see an article about estimating the average,
variance, skewness, and kurtosis, for the number of occurrences (as
consecutive submatrix) of the matrix [[1, 0, 2], [1, 0, 2], [1, 0, 1]], in
row-echelon matrices of dimension $k\times 2k$ over $GF(3)$ for $k$ from $50$
to $60$ by simulating $1000$ times, look at the output file

http://www.math.rutgers.edu/~zeilberg/tokhniot/oCalabiWilf1a.txt .

$\bullet$ If you want to see an article about estimating the average,
variance, skewness, and kurtosis of the number of $1$s in $k\times 2k$ row-
echelon matrices over $GF(3)$ for $50\leq k\leq 55$ , by simulating $1000$
times (and repeating each run three times, to compare notes), and comparing
the estimated values to the (theoretically computed) exact values, look at

http://www.math.rutgers.edu/~zeilberg/tokhniot/oCalabiWilf2.txt .

$\bullet$ If you want to see an article about estimating the average of the
number of $1$s in $k\times 3k$ row-echelon matrices over $GF(2)$ for $100\leq
k\leq 110$ , by simulating $1000$ times (and repeating each run three times,
to compare notes), and comparing the estimated average to the (theoretically
computed) exact value, look at

http://www.math.rutgers.edu/~zeilberg/tokhniot/oCalabiWilf3.txt .

$\bullet$ If you want to see an article on the average minimal weight of
vectors in $k$-dimensional subspaces of $GF(2)^{n}$ for $k$ from $1$ to, $5$,
and $n$ from $10$ to $100$, in increments of $10$, look at the output file

http://www.math.rutgers.edu/~zeilberg/tokhniot/oCalabiWilf4.txt .

For more output files see the web-page of this paper:

https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/calabi.html .

Readers are welcome to experiment with our Maple package to their heart’s
content, and generate many more output files.

Conclusion: Rest in peace, my dear heroes, Herb and Gene, you both did great
work on your own, and in collaboration with other people, in very different
parts of mathematics, but this sole joint work is also great!

References

[CW] Eugenio Calabi and Herbert S. Wilf, On the Sequential and Random
Selection of Subspaces over a Finite Field, J. Of Combinatorial Theory (Ser.
A) 22 (1977), 107-109. Reprinted in: Eugenio Calabi, “Collected works”, Edited
by Jean-Pierre Bourguignon, Xiuxiong Chen and Simon Donaldson. With
contributions by Shing-Tung Yau, Blaine Lawson, Marcel Berger and Claude
LeBrun Springer, Berlin, (2020), 481-483.
https://sites.math.rutgers.edu/~zeilberg/akherim/CalabiWilf1977.pdf .

[NW] Albert Nijenhuis and Herbert S. Wilf, “Combinatorial Algorithms”,
Academic Press. First edition: 1975. Second edition: 1978.
https://www2.math.upenn.edu/~wilf/website/CombinatorialAlgorithms.pdf .

[P] Igor Pak, When and how n choose k, Randomization methods in algorithm
design, 191-238, 1997. https://www.math.ucla.edu/~pak/papers/nk13.pdf .

[W] Herbert S. Wilf, A unified setting for sequencing, ranking, and selection
algorithms for combinatorial objects, Advances in Math 24 (1977) , 281-291.
https://www2.math.upenn.edu/~wilf/website/Unified%20setting.pdf .

Shalosh B. Ekhad, c/o D. Zeilberger, Department of Mathematics, Rutgers
University (New Brunswick), Hill Center-Busch Campus, 110 Frelinghuysen Rd.,
Piscataway, NJ 08854-8019, USA. Email: ShaloshBEkhad at gmail dot com .

Doron Zeilberger, Department of Mathematics, Rutgers University (New
Brunswick), Hill Center-Busch Campus, 110 Frelinghuysen Rd., Piscataway, NJ
08854-8019, USA. Email: DoronZeil at gmail dot com .

Oct. 27, 2023.


# An Energy-driven Network Function Virtualization for Multi-domain Software
Defined Networks

Kuljeet Kaur1, Member, IEEE, Sahil Garg1, Member, IEEE, Georges Kaddoum1,
Member, IEEE,
François Gagnon1, Senior Member, IEEE, Neeraj Kumar2, Senior Member, IEEE, and
Syed Hassan Ahmed3, Senior Member, IEEE
(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>and<EMAIL_ADDRESS> 1Electrical Engineering Department, École de
technologie supérieure, Montréal, QC H3C 1K3, Canada.
2CSED, Thapar Institute of Engineering & Technology, Patiala (Punjab), India.
3Department of Computer Science, Georgia Southern University, Statesboro, GA
30460, USA.

###### Abstract

Network Functions Virtualization (NFV) in Software Defined Networks (SDN)
emerged as a new technology for creating virtual instances for smooth
execution of multiple applications. Their amalgamation provides flexible and
programmable platforms to utilize the network resources for providing Quality
of Service (QoS) to various applications. In SDN-enabled NFV setups, the
underlying network services can be viewed as a series of virtual network
functions (VNFs) and their optimal deployment on physical/virtual nodes is
considered a challenging task to perform. However, SDNs have evolved from
single-domain to multi-domain setups in the recent era. Thus, the complexity
of the underlying VNF deployment problem in multi-domain setups has increased
manifold. Moreover, the energy utilization aspect is relatively unexplored
with respect to an optimal mapping of VNFs across multiple SDN domains. Hence,
in this work, the VNF deployment problem in multi-domain SDN setup has been
addressed with a primary emphasis on reducing the overall energy consumption
for deploying the maximum number of VNFs with guaranteed QoS. The problem in
hand is initially formulated as a “Multi-objective Optimization Problem” based
on Integer Linear Programming (ILP) to obtain an optimal solution. However,
the formulated ILP becomes complex to solve with an increasing number of
decision variables and constraints with an increase in the size of the
network. Thus, we leverage the benefits of the popular evolutionary
optimization algorithms to solve the problem under consideration. In order to
deduce the most appropriate evolutionary optimization algorithm to solve the
considered problem, it is subjected to different variants of evolutionary
algorithms on the widely used MOEA framework (an open source java framework
based on multi-objective evolutionary algorithms). The experimental results
demonstrate that the proposed scheme achieves better results in comparison to
the $\epsilon$-Non-dominated Sorting Genetic Algorithm (NSGA)-II
($\epsilon$-NSGA-II) with the respect to the overall energy consumption and
optimal deployment of VNFs in multi-domain SDN scenarios.

###### Index Terms:

Network Function Virtualization, Software Defined Network, Evolutionary
Optimization, Energy Consumption, Multi-objective Optimization

##  I Introduction

With the development of new agile technologies such as 5G (and beyond), the
number of connected devices are increasing at a rapid pace. However, the lack
of direct application program interface (API) interactions with network
services impose limitations on the handling of the rapidly changing traffic
patterns. Also, the rise of cloud, Internet of Things (IoT) and AI-driven
applications have made this task more challenging. To cope with the ever-
increasing demands of networks in a cost-efficient manner, the conventional
networks are demanding more flexible and dynamic solutions [1, 2]. Towards
this end, network softwarization is gaining huge popularity which supports
rapid provisioning and deployment of network services to meet QoS requirements
[3]. By reducing overhead, automating infrastructure and increasing
maintainability, it aims to build agile and flexible networks to provision
convergent networking solutions.
In this direction, Software Defined Networks (SDN) and Network Functions
Virtualization (NFV) have been introduced [4]. SDN is an exciting and evolving
paradigm that turns the network into a flexible and programmable platform able
to optimally utilize the resources. It separates the control plane of a
network from the data plane to make it centrally manageable [5]. NFV, on the
other hand, represents the network service as a series of virtual network
functions (VNFs). In NFV architectures, individual VNFs can be chained
together to perform the desired sequence, known as service function chain
(SFC). Using SFC, a large number of VNFs can be connected together in an NFV
environment. To support a fully virtualized network, VNFs provide network
scalability and agility, while enabling better provisioning of network
resources [6].
To meet the dynamic requirements of the network resources, several solutions
have been proposed in the literature for SFCs deployment. For example, Chen et
al. [7] proposed a fully decentralized online approach for service chaining
based on stochastic dual-gradient method. Similarly, Sun et al. [8] studied
the problem of SFC orchestration using a mesh aggregation approach where a
feedback mechanism was deployed to improve its success rate. In another work,
Sun et al. [9] proposed a heuristic algorithm using a restricted Boltzmann
machine and cross-entropy approach for low-latency applications in NFV
networks. Likewise, Gupta et al. [10] proposed a VNF service chain placement
model to facilitate the dynamic service demands of network traffic. Despite
the rapid progression, the dynamic composition of services is still facing
imminent challenges which need immediate attention. In this direction, the
integration of SDN and NFV has been proposed as a viable solution; where NFV
provides the basic functionality to the network, while SDN controls and
orchestrates them for specific uses. Thus, the convergence of both
technologies facilitates the flexible deployment of networks and service
delivery over them.
In the similar context, Medhat et al. [11] explored the limitations of current
SFC approaches in next-generation networks and showed that approaches
integrating SDN and NFV technologies have higher SFC scalability and
flexibility compared to others. To ensure the benefits of the service-chaining
policy, Ding et al. [12] proposed OpenSCaaS platform by integrating SDN and
NFV. Furthermore, Zhang et al. [13] presented a framework to construct
efficient SFCs by combining SDN and NFV. Bruschi et al. [14] proposed a
softwarized SDN/NFV approach to support the cloud services in a more scalable
and sustainable way. Likewise, Pie et al. [15] developed a novel routing
algorithm to solve the SFC problem in SDN and NFV-enabled network. In order to
minimize the resource consumption flow costs with SFC requests, binary integer
programming was used. Although several attempts have been made to tackle the
service deployment of networks, they all rely on single domain SDN networks.
However, the rapidly growing network traffic and constant enlargement of
network scales pose challenges to single domain SDN architectures in terms of
latency, overhead, fault tolerance and load balancing, thereby affecting the
overall QoS experience for end-users.
To overcome these challenges, large-scale SDN networks have been proposed.
However, the deployment of network services on a large-scale SDN network
involves the distribution of VNFs on multiple SDN domains. Thus, the service
provisioning across these networks requires the collaboration among different
SDN controllers which results in additional network overhead. However, the
deployment of network services across multi-domain SDNs is still in its early
phases [16]. Further, the deployment of network services in multi-domain SDN
networks requires an efficient allocation of resources for different service
requests at a reduced cost. Thus, it is of primary interest to implement
energy efficient solutions which dynamically adapt the usage of resources
according to the demands.

###  I-A Contributions

This paper proposes an energy-efficient approach to deploy NFVs in multi-
domain SDN networks. To the best of our knowledge, the evaluation of energy
characteristics in multi-domain SDN setups is a novel research initiative and
has not been explored. The key contributions of the proposed wok are
summarized as follows:

  1. 1.

We analyze the problem of VNF deployment in multi-domain SDN environments
(associated with energy consumption and resource allocation challenges), and
mathematically formulate a multi-objective optimization problem for the
effective deployment of NFV in SDN. The formulated problem aims at achieving
an optimal mapping of VNFs across a multi-domain SDN setup while minimizing
the overall energy consumption level.

  2. 2.

We then evaluate the proposed multi-objective problem on different
evolutionary algorithms to have an optimal solution in polynomial time using
the widely used MOEA framework. To attain this objective, different variants
of Non-dominated Sorting Genetic Algorithm (NSGA) are used namely NSGA-II,
NSGA-III and $\epsilon-$NSGA-II.

  3. 3.

We then extensively assess the performance of the proposed scheme using the
most appropriate evolutionary algorithm under different scenarios in terms of
number of VNFs, number of SDN domains, number of constraints and decision
variables.

###  I-B Organization

The rest of this paper is organized as follows. Section II presents the system
model for the proposed scheme. In Section III, the problem of VNF deployment
in multi-domain SDN is formulated. Finally, in Section IV, the formulated
problem is subjected to different variants of NSGA in different networks
scenarios. Finally, this work is concluded in Section V with future
directions.

##  II System Model

In the following section, we detail the system model for the next-generation
communication systems empowered by SDN and NFV. The designed solution helps in
the optimized mapping of the VNFs to the available SDN nodes (scattered across
multiple domains) keeping the energy consumption and resource utilization
profile in check.

Figure 1: Architectural model based on SDN and NFV.

Fig. 1 depicts the system model of the proposed communication system. It
represents the coupling of SDN with NFV for seamless networking operations.
SDN is a software-based programming paradigm for modern networks that rely on
software reconfigurations for changing the network policies. As shown in the
figure, the SDN platform is comprised of a data plane and a control plane.
Here, the data plane consists of the Open-Flow (OF)-enabled routers and
switches with dedicated flow table entries. These OF-flow tables are primarily
controlled by the SDN’s control plane via the southbound API, i.e., OpenFlow
protocols. In other words, the policy changes are assumed by the SDN’s
controller; which are reproduced at the data plane [17].

On the other hand, NFV is an upcoming networking paradigm that helps to
segregate the network functions from the underlying proprietary hardware.
Further, it helps to execute the network functions on software with flexible
implementation which are referred to as VNFs. By coupling NFV with the
powerful paradigm of SDN, the network services can be established with ease by
redirecting the traffic through the appropriate VNFs in a sequential manner.
The Internet Engineering Task Force (IETF) defines the abstract notion of VNFs
using SFC–a means to implement VNFs. These SFCs are managed by the SDN’s
centralized controller and deployed on the underlying hardware. With the SDN-
NFV coupling, the deployment of SFCs becomes flexible without the need to
change the hardware.

Figure 2: Deployment of VNFs across SDN domains.

However, the constant enlargement of networks has resulted in wide-scale SDNs
with several interconnected domains referred to as multi-domain SDNs [16]. All
the domains are individually administrated by their respective controllers and
are linked with each other to provide the required network services. The
provision of these network services involves the successful deployment of VNFs
(belonging to different SCFs) across different SDN domain nodes (physical or
virtual). Nevertheless, the collaboration between different SDN domain
controllers is essential for the efficient distribution of resources amongst
the VNFs. The same can be achieved using east and west-bound APIs. In summary,
the focus of the proposed work is to manage the underlying resources
effectively while attaining minimal energy consumption. In the proposed work,
the same is achieved by mapping VNFs to the available array of SDN domains
while assigning the required resources (CPU cores, memory and storage) to VNFs
for their execution. The details of the overall deployment process are
illustrated using Fig. 2.

##  III Problem Formulation

In this section, the problem of VNF deployment across a multi-domain SDN set-
up is illustrated. The key objective functions of the considered problem
followed by their corresponding constraints are presented.

A. List of Objective Functions

As discussed earlier, the proposed scheme supports two major objective
functions ($\mathfrak{F}(\mathbb{X}_{ij})$) as illustrated under. Here,
$\mathbb{X}_{ij}$ denotes the decision variable expressed as follows:

$\displaystyle\mathbb{X}_{ij}=$ $\displaystyle\begin{cases}1:\text{If $i^{th}$
VNF is deployed on the $j^{th}$ SDN domain}\\\ 0:\text{Otherwise}\end{cases}$
(1)

In the above equation, the nature of $\mathbb{X}_{ij}$ is binary, i.e.,
$\mathbb{X}_{ij}\in\\{0,1\\}$. Here, $\mathbb{X}_{ij}=1$ indicates that the
$i^{th}$ VNF is deployed successfully on the $j^{th}$ SDN domain, while
$\mathbb{X}_{ij}=0$ suggests that $i^{th}$ VNF is not deployed on the $j^{th}$
SDN domain.

1) Maximum deployment of VNFs: One of the foremost objectives of the proposed
scheme is to ensure maximum deployment of VNFs on the available array of SDN
nodes. This is attained by using a maximization operation over the summation
of $\mathbb{X}_{ij}$ as follows:

$\mathfrak{F}_{1}(\mathbb{X}_{ij})=\max\sum_{i=1}^{\mathrm{M}}\sum_{j=1}^{\mathrm{N}}\mathbb{X}_{ij}$
(2)

In the above equation, $\mathrm{M}$ denotes the total number of VNFs to be
deployed.

2) Minimum Energy Consumption: This objective function deals with the
minimization of the overall energy consumption involved in executing VNFs on
the available SDN domains. It is denoted as follows:

$\mathfrak{F}_{2}(\mathbb{X}_{ij})=\min\sum_{i=1}^{\mathrm{M}}\sum_{j=1}^{\mathrm{N}}\mathbb{X}_{ij}\times\mathbb{E}_{ij}$
(3)

where, $\mathrm{N}$ denotes the total number of SDN domains and variable
$\mathbb{E}_{ij}$ refers to the energy consumption (in kWh) associated with
the execution of the $i^{th}$ VNF over the $j^{th}$ SDN domain.

B. List of Constraints

This section presents different set of constraints that should be respected
while achieving energy-aware deployment of VNFs in the multi-domain SDN setup.

1) VNF deployment constraint: The following constraint ensures that only a
single VNF is allocated at most a single SDN domain as expressed using the
below mentioned equation. The said constraint restricts the multiple
deployment of an individual VNF in multi-domain SDN environments and helps to
minimize the overall energy consumption.

$\displaystyle\mathrm{C}_{1}:$
$\displaystyle\mathbb{X}_{i1}+\mathbb{X}_{i2}+\dots+\mathbb{X}_{i\mathrm{N}}=1;$
$\displaystyle\forall i,$ $\displaystyle\implies$
$\displaystyle\sum_{j=1}^{\mathrm{N}}\mathbb{X}_{ij}=1;$ $\displaystyle\forall
i.$ (4)

2) Resource restriction with respect to CPU cores: The following constraint
ensures that VNF deployment on different physical nodes always meets their
required CPU core needs. The said restriction helps to allocate the requested
number of cores to the mapped VNFs using the following equation. In other
words, the constraints ensure that no more VNFs can be mapped to a physical
node without the required number of CPU cores requested by the corresponding
VNF.

$\displaystyle\mathrm{C}_{2}:$
$\displaystyle\sum_{i=1}\mathbb{X}_{1j}\times\mathbb{C}_{1}\dots+\mathbb{X}_{\mathrm{M}j}\times\mathbb{C}_{\mathrm{M}}\leq\mathfrak{C}_{j};$
$\displaystyle\forall j,$ $\displaystyle\implies$
$\displaystyle\sum_{i=1}^{\mathrm{M}}\mathbb{X}_{ij}\times\mathbb{C}_{i}\leq\mathfrak{C}_{j};$
$\displaystyle\forall j.$ (5)

In the above equation, $\mathbb{C}_{i}$ refers to the number of the CPU cores
required by the $i^{th}$ VNF and $\mathfrak{C}_{j}$ denotes the upper limit on
the $j^{th}$ node’s CPU core capacity.

3) Resource restriction with respect to Memory: In line, with the above
mentioned constraint, the following restriction imposes the upper limit of the
main memory that could be allocated to the VNFs on a particular SDN domain
node.

$\displaystyle\mathrm{C}_{3}:$
$\displaystyle\sum_{i=1}\mathbb{X}_{1j}\times\mathbb{M}_{1}\dots+\mathbb{X}_{\mathrm{M}j}\times\mathbb{M}_{\mathrm{M}}\leq\mathfrak{M}_{j};$
$\displaystyle\forall j,$ $\displaystyle\implies$
$\displaystyle\sum_{i=1}^{\mathrm{M}}\mathbb{X}_{ij}\times\mathbb{M}_{i}\leq\mathfrak{M}_{j};$
$\displaystyle\forall j.$ (6)

here, variables $\mathbb{M}_{i}$ and $\mathfrak{M}_{j}$ represent the main
memory required by the $i^{th}$ VNF and the upper limit on the main memory
capacity of the $j^{th}$ node, respectively.

4) Resource restriction with respect to Storage: The following equation
represents the resource restriction corresponding to disk storage allocation.

$\displaystyle\mathrm{C}_{4}:$
$\displaystyle\sum_{i=1}\mathbb{X}_{1j}\times\mathbb{S}_{1}\dots+\mathbb{X}_{\mathrm{M}j}\times\mathbb{S}_{\mathrm{M}}\leq\mathfrak{S}_{j};$
$\displaystyle\forall j,$ $\displaystyle\implies$
$\displaystyle\sum_{i=1}^{\mathrm{M}}\mathbb{X}_{ij}\times\mathbb{S}_{i}\leq\mathfrak{S}_{j};$
$\displaystyle\forall j.$ (7)

The variable $\mathbb{S}_{i}$ defined in the above equation denotes the disk
storage requested by the $i^{th}$ VNF; while the variable $\mathfrak{S}_{j};$
denotes the upper limit on disk storage for the $j^{th}$ node.

C. Overall Problem at hand

The overall problem can be summarized as follows:

$\displaystyle\mathfrak{F}(\mathbb{X}_{ij})=$
$\displaystyle\big{[}-\mathfrak{F}_{1}(\mathbb{X}_{ij}),$
$\displaystyle\mathfrak{F}_{2}(\mathbb{X}_{ij})\big{]}$ (8) such that
Constraints in
$\displaystyle\begin{cases}\text{~{}Eq.~{}}\eqref{eq:Constraint1}\\\
\text{~{}Eq.~{}}\eqref{eq:Constraint4}\\\
\text{~{}Eq.~{}}\eqref{eq:Constraint2}\\\
\text{~{}Eq.~{}}\eqref{eq:Constraint3}\\\ \end{cases}$ holds

As detailed above, the considered problem comprises of the two objective
functions $\mathfrak{F}_{1}(\mathbb{X}_{ij})$ and
$\mathfrak{F}_{2}(\mathbb{X}_{ij})$. The minus sign before
$\mathfrak{F}_{1}(\mathbb{X}_{ij})$ denotes maximization operation, i.e., the
aim is to maximize the deployment of VNFs. On the other hand,
$\mathfrak{F}_{2}(\mathbb{X}_{ij})$ depicts the minimization operation over
the energy consumption objective function. The associated constraints span
from Eq. (4) till Eq. (7). It is worth noticing here that the nature of the
problem mentioned above is linear with linear constraints; which suggests that
the problem under consideration is ILP with binary decision variable
$\mathbb{X}_{ij}$. However, solving an ILP with a huge number of parameters
(say thousands of $\mathbb{X}_{ij}$) is quite a gigantic task. This is because
with increasing Internet traffic, multi-domains SDNs are expected to grow
exponentially, which in turn would lead to a substantial increase in the
number of VNFs ($\mathrm{M}$) and SDN-domains ($\mathrm{N}$). Consider the
following example,

$\displaystyle\mathrm{M}$ $\displaystyle=100;$ $\displaystyle\mathrm{N}$
$\displaystyle=20;$ $\displaystyle\therefore\text{~{}No. of Decision
Variables}$ $\displaystyle~{}(|\mathbb{X}_{ij}|)=\mathrm{M}\times\mathrm{N}$
$\displaystyle=100*20=2000$ No. of Constraints
$\displaystyle~{}(|\mathrm{C}_{k}|)=\mathrm{M}+3\times\mathrm{N}$
$\displaystyle=100+3*20=160$

Thus, with the constant increase in parameters, an optimal solution for the
above problem would be time-consuming and impossible in polynomial time.
Henceforth, the proposed scheme leverages the advantages of a popular multi-
objective evolutionary algorithm called $\epsilon$-NSGA-II. It is
characterized by its fast convergence, low overall complexity, explicit
diversity preservation capability, epsilon-dominance, adaptive population
sizing, and automatic termination [18, 19].

The experimental evaluation against the existing state-of-the-art multi-
objective evolutionary algorithms demonstrate that $\epsilon$-NSGA-II is best
suited for our problem, as justified by the preliminary results demonstrated
in the upcoming section.

##  IV Results and Discussions

This section illustrates the experimental evaluation of the proposed scheme
compared to different state-of-the-art schemes on different evaluation
parameters.

To evaluate the proposed VNF deployment, the widely accepted MOEA Framework
has been employed [20]. It is an open source Java library which is heavily
used by the research community for designing and evaluating different multi-
objective evolutionary algorithms such as $\epsilon$-NSGA-II, NSGA-II, NSGA-
III, generalized differential evolution (GDE3), strength pareto evolutionary
algorithm 2 (SPEA2), indicator-based evolutionary algorithm (IBEA), speed-
constrained multi-objective particle swarm optimization (SMPSO) covariance
matrix adaptation evolution strategy (CMA-ES), multi-objective evolutionary
algorithm with decomposition (MOEA/D), and $\epsilon$-MOEA. Additionally, it
also supports rapid design and development of customized problems and
algorithms. For the proposed scheme, the formulated multi-objective
optimization problem (defined Eq. (III)) has been scripted in java using MOEA
Framework and has been evaluated with existing state-of-the-art algorithms
such as NSGA-II, NSGA-III, and $\epsilon$-NSGA-II.

(a) Pareto front evaluation.

(b) Elapsed time evaluation.

Figure 3: An illustrative comparison of NSGA-II, NSGA-III and $\epsilon$-NSGA-
II for solving the problem under consideration.

The experiment evaluation for the proposed multi-objective scheme was carried
out across two phases (referred to as Case Studies). In the first phase of the
assessment, the proposed scheme was subjected to different multi-objective
evolutionary algorithms, and their relative performance was compared
considering Pareto front and elapsed time. Based on the comparison, it was
concluded that $\epsilon$-NSGA-II was the most appropriate algorithm to
address the proposed NFV deployment across multi-domain SDN setups.
Accordingly, in the next phases, the efficacy of the proposed scheme was
evaluated using $\epsilon$-NSGA-II on a dataset spanning different scenarios.
The detailed description of these phases is illustrated below.

(a) Number of VNFs to be deployed.

(b) Number of SDN domain nodes.

(c) Execution of $\mathfrak{F}_{1}$ for successful deployment of VNFs.

(d) Execution of $\mathfrak{F}_{2}$ for minimizing energy consumption.

Figure 4: A illustration of executing the proposed multi-objective
optimization problem for VNF deployment in multi-domain SDN using
$\epsilon$-NSGA-II.

A. Case Study I

The problem of VNF deployment in multi-domain environments has been evaluated
using different versions of NSGA namely NSGA-II, NSGA-III, and
$\epsilon$-NSGA-II. The related results are summarized in Fig. 3. Each of
these algorithms has been evaluated with 10,000 function evaluations with 20
seeds in the MOEA framework’s Diagnostic Tool. The related results regarding
pareto front and elapsed time are summarized in Fig. 3(a) and Fig. 3(b),
respectively. The pareto front evaluations show that $\epsilon$-NSGA-II
achieves the most optimal performance followed by NSGA-II and NSGA-III; in
terms of trade-off between the two objective function
$\mathfrak{F}_{1}(\mathbb{X}_{ij})$ and $\mathfrak{F}_{2}(\mathbb{X}_{ij})$.
Although the performance of $\epsilon$-NSGA-II and NSGA-II is almost
identical; however the former results in much less energy consumption
deployments as compared to the other. Additionally, the computational time of
the three algorithms is summarized in Fig. 3(b). Similar to the previous case,
$\epsilon$-NSGA-II achieves the optimal performance with a reduced value of
elapsed time. This suggests that $\epsilon$-NSGA-II achieves the optimal
performance with quick computational processing.

Thus, for both parameters, $\epsilon$-NSGA-II stands out which makes it the
ideal choice to solve the considered VNF deployment challenge in a multi-
domain SDN setup.

B. Case Study II

During this case study, the performance of the proposed multi-objective
approach for VNF deployment has been evaluated using $\epsilon$-NSGA-II under
different scenarios in terms of number of VNFs deployed, number of SDN domain
nodes, number of decision variables and number of constraints. The variability
concerning the first two parameters can be understood from the data plotted in
Figs. 4(a) and Fig 4(b). For the other two parameters, the following equations
can be used for the computational purpose.

No. of Decision Variables $\displaystyle~{}(|\mathbb{X}_{ij}|)=$
$\displaystyle\mathrm{M}\times\mathrm{N}$ No. of Constraints
$\displaystyle~{}(|\mathrm{C}_{k}|)=$
$\displaystyle\mathrm{M}+3\times\mathrm{N}$

In order to achieve successful deployment of VNFs in a multi-domain setup, Eq.
(8) was evaluated using $\epsilon$-NSGA-II. The related results in terms of
successful VNF deployment and energy consumptions statistics are highlighted
in Fig. 4(c) and Fig. 4(d), respectively. The obtained results evidently
support the fact that the proposed scheme achieves an optimal trade-off
between the two competing objective functions
$\mathfrak{F}_{1}(\mathbb{X}_{ij})$ and $\mathfrak{F}_{2}(\mathbb{X}_{ij})$.
Overall, in comparison with the existing NSGA variants, the proposed approach
using $\epsilon$-NSGA-II gives better performance regarding energy consumption
and successful VNF deployment.

##  V Conclusion

Gradual evolutions in the networking infrastructures have resulted in the
convergence of SDN and NFV to provide seamless deployment of network services.
Additionally, this convergence has led to the amplification of SDNs from
single domain to multiple domains. Under such multi-domain SDNs, the
realization of NFV would become a challenging task particularly due to the
complex interactions between SDN controllers of different domains.
Particularly, the energy consumption management will also evolve as a
challenging task. Thus, in this paper, the problem of NFV deployment in multi-
domain SDNs has been addressed. The main focus of the proposed work was to
minimize the energy consumption caused by VNF deployment across multiple
domain SDN while achieving their optimal deployment with the required number
of resources. In summary, the said problem has been modeled as an ILP and
solved using different variants of NSGA-a popular evolutionary optimization
algorithm. However, experimental evaluations suggests that $\epsilon$-NSGA-II
gives the optimal trade-off between the considered objective functions and is
best suited for the problem under consideration. Additionally, the
experimental results also suggest that the proposed VNF deployment solution
also reduces the energy consumption to a large extent.

In the near future, we will extend this work to locality-based VNF deployment;
where VNFs of a particular SFC would be deployed in a specific SDN domain.

## References

  * [1] A. Singh, S. Garg, S. Batra, N. Kumar, and J. J. Rodrigues, “Bloom filter based optimization scheme for massive data handling in IoT environment,” _Future Generation Computer Systems_ , vol. 82, pp. 440–449, 2018.
  * [2] S. Garg, A. Singh, K. Kaur, S. Batra, N. Kumar, and M. S. Obaidat, “Edge-Based Content Delivery for Providing QoE in Wireless Networks Using Quotient Filter,” in _IEEE International Conference on Communications (ICC), Kansas City, USA_ , May 2018.
  * [3] K. Kaur, S. Garg, G. S. Aujla, N. Kumar, J. J. Rodrigues, and M. Guizani, “Edge Computing in the Industrial Internet of Things Environment: Software-Defined-Networks-Based Edge-Cloud Interplay,” _IEEE Communications Magazine_ , vol. 56, no. 2, pp. 44–51, 2018.
  * [4] V. Nguyen, A. Brunstrom, K. Grinnemo, and J. Taheri, “SDN/NFV-Based Mobile Packet Core Network Architectures: A Survey,” _IEEE Communications Surveys Tutorials_ , vol. 19, no. 3, pp. 1567–1602, 2017.
  * [5] S. Garg, K. Kaur, N. Kumar, and J. J. P. C. Rodrigues, “Hybrid Deep-Learning-Based Anomaly Detection Scheme for Suspicious Flow Detection in SDN: A Social Multimedia Perspective,” _IEEE Transactions on Multimedia_ , vol. 21, no. 3, pp. 566–578, 2019.
  * [6] R. Mijumbi, J. Serrat, J. Gorricho, N. Bouten, F. D. Turck, and R. Boutaba, “Network Function Virtualization: State-of-the-Art and Research Challenges,” _IEEE Communications Surveys Tutorials_ , vol. 18, no. 1, pp. 236–262, 2016.
  * [7] X. Chen, W. Ni, T. Chen, I. Collings, X. Wang, R. P. Liu, and G. B. Giannakis, “Multi-Timescale Online Optimization of Network Function Virtualization for Service Chaining,” _IEEE Transactions on Mobile Computing_ , 2018, DOI: 10.1109/TMC.2018.2885301.
  * [8] G. Sun, Y. Li, D. Liao, and V. Chang, “Service Function Chain Orchestration Across Multiple Domains: A Full Mesh Aggregation Approach,” _IEEE Transactions on Network and Service Management_ , vol. 15, no. 3, pp. 1175–1191, 2018.
  * [9] G. Sun, G. Zhu, D. Liao, H. Yu, X. Du, and M. Guizani, “Cost-Efficient Service Function Chain Orchestration for Low-Latency Applications in NFV Networks,” _IEEE Systems Journal_ , pp. 1–13, 2018, DOI: 10.1109/JSYST.2018.2879883.
  * [10] A. Gupta, M. F. Habib, U. Mandal, P. Chowdhury, M. Tornatore, and B. Mukherjee, “On Service-Chaining Strategies using Virtual Network Functions in Operator Networks,” _Computer Networks_ , vol. 133, pp. 1–16, 2018.
  * [11] A. M. Medhat, T. Taleb, A. Elmangoush, G. A. Carella, S. Covaci, and T. Magedanz, “Service function chaining in next generation networks: State of the art and research challenges,” _IEEE Communications Magazine_ , vol. 55, no. 2, pp. 216–223, 2017.
  * [12] W. Ding, W. Qi, J. Wang, and B. Chen, “OpenSCaaS: an open service chain as a service platform toward the integration of SDN and NFV,” _IEEE Network_ , vol. 29, no. 3, pp. 30–35, 2015.
  * [13] J. Zhang, Z. Wang, N. Ma, T. Huang, and Y. Liu, “Enabling Efficient Service Function Chaining by Integrating NFV and SDN: Architecture, Challenges and Opportunities,” _IEEE Network_ , vol. 32, no. 6, pp. 152–159, 2018.
  * [14] R. Bruschi, F. Davoli, P. Lago, A. Lombardo, C. Lombardo, C. Rametta, and G. Schembra, “An SDN/NFV Platform for Personal Cloud Services,” _IEEE Transactions on Network and Service Management_ , vol. 14, no. 4, pp. 1143–1156, 2017.
  * [15] J. Pei, P. Hong, K. Xue, and D. Li, “Resource Aware Routing for Service Function Chains in SDN and NFV-Enabled Network,” _IEEE Transactions on Services Computing_ , 2018, DOI: 10.1109/TSC.2018.2849712.
  * [16] C. Zhang, X. Wang, F. Li, M. Huang, and Q. He, “Network service chains deployment across multiple SDN domains,” _International Journal of Communication Systems_ , vol. 31, no. 18, p. e3826, 2018.
  * [17] R. Chaudhary, N. Kumar, and S. Zeadally, “Network service chaining in fog and cloud computing for the 5G environment: Data management and security challenges,” _IEEE Communications Magazine_ , vol. 55, no. 11, pp. 114–122, 2017.
  * [18] P. Reed, B. S. Minsker, and D. E. Goldberg, “Simplifying multiobjective optimization: An automated design methodology for the nondominated sorted genetic algorithm-II,” _Water Resources Research_ , vol. 39, no. 7, 2003\.
  * [19] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: NSGA-II,” _IEEE transactions on evolutionary computation_ , vol. 6, no. 2, pp. 182–197, 2002.
  * [20] D. Hadka, “MOEA Framework: A Free and Open Source Java Framework for Multiobjective Optimization,” [Accessed on: Nov. 2018]. [Online]. Available: http://moeaframework.org


# Note on stability and holographic subregion complexity

Mohammad<EMAIL_ADDRESS>Department of Physics, Shahid
Beheshti University, 1983969411, Tehran, Iran  Mahsa<EMAIL_ADDRESS>Department of Physics, Shahid Beheshti University, 1983969411, Tehran, Iran

###### Abstract

We study holographic subregion complexity in a spatially anisotropic field
theory, which expresses a confinement-deconfinement phase transition. Its
holographic dual is a five-dimensional anisotropic holographic model
characterized by a Van der Waals-like phase transition between small and large
black holes. We propose a new interpretation from the informational
perspective to determine the stable and unstable thermodynamically solutions.
According to this proposal, the states which need (more) less information to
be specified characterize the (un) stable solutions. We similarly offer an
interpretation to determine the stable and unstable solutions based on the
resource of a computational machine, such that the solutions are (un) stable
if computational resource (decreases) increases with the increase of
temperature. We observe that the effect of anisotropy on holographic subregion
complexity is decreasing. This decreasing effect can be interpreted by
considering a whole closed system consisting of the state and its environment
in which the complexity of the mixed state decreases and complexity of the
environment increases.

###### Contents

  1. I Introduction
  2. II Review on the background
  3. III Holographic subregion complexity
  4. IV Numerical results
    1. IV.1 Isotropic background
    2. IV.2 Anisotropic background
      1. IV.2.1 Parallel case
      2. IV.2.2 Perpendicular case

##  I Introduction

Gauge-gravity duality or more generally the holographic idea, is an
alternative and useful tool for investigation of a strong coupling system,
where standard methods are powerless CasalderreySolana:2011us . Quark-gluon
plasma (QGP) produced in heavy ion collisions presents a strongly coupled
fluid with a small viscosity that is an example of such a system one . It is
believed that the QGP is produced after a very short time after the collision,
$\tau_{\rm{therm}}\approx\rm{few}~{}0.1~{}\rm{fm/c}$ and there are indications
that the anisotropic stage of the QGP takes place for time scales of
$0.1~{}\rm{fm/c}\mathchar 13358\relax\tau\mathchar 13358\relax
0.3-2~{}\rm{fm/c}$ two .

According to gauge-gravity duality, a strongly coupled gauge theory in
$d$-dimensional space-time corresponds to a classical gravity in
$d+1$-dimensional space-time, therefore, to describe different phenomena in
strongly coupled and anisotropic field theories in this case, it suffices to
find a proper description on the gravity side. Different effects of anisotropy
have been studied in holographic contexts three ; four . The Lifshitz-like
background, which is holographically dual to an anisotropic plasma, is applied
to study various aspects of the anisotropic plasma. For example, in these
models, holographic estimates of the total multiplicity for certain values of
critical exponent make what is desirable from experimental data five .

One of the interesting areas of this duality is a great connection which is
developed to connect the quantities in quantum information theory and certain
geometric quantities in the bulk geometry. Complexity is an important concept
in quantum information theory. This quantity refers to the time and space
resources required to perform a computation efficiently co1 . In the context
of quantum field theory, complexity is defined as the minimum number of
unitary operators needed to prepare a target state from a reference state co2
. In order to describe the complexity holographically, the complexity of pure
states in the whole boundary field theory has two proposals namely CV
(complexity=volume) conjecture and CA (complexity=action) conjecture. In CA
conjecture, the complexity is obtained by the bulk action evaluated on the
Wheeler-de Witt patch anchored at some boundary time. In CV conjecture, the
complexity is defined as the volume of the extremal/maximal volume of a
codimensional-one hypersurface in the bulk ending on a time slice of the
boundary co3 ; co4 . Inspired by the Hubney-Ryu-Takayanagi proposal, the CV
proposal extends to be defined on subregions corresponding to the complexity
of mixed states, which is known as holographic subregion complexity (HSC), in
which the complexity of a subsystem on the boundary equals the volume of
codimensional-one hypersurface enclosed by Hubney-Ryu-Takayanagi surface
alishahiha . There are a lot of works on CV, CA and HSC for various gravity
models in the literature co5 ; co6 ; co7 ; co8 ; co9 ; co10 ; co11 ; co12 ;
co13 ; co14 ; co15 .

According to thermodynamics of complexity, complexity is like the entropy of
an auxiliary classical system, which describes the evolution of unitary
operators unco . For instance, the second law of complexity is the usual
second law of thermodynamics applied to the auxiliary classical system.
According to this analogy, the difference between the maximum entropy and the
actual entropy, which is a resource for doing a work has a parallel definition
in the context of complexity known as uncomplexity. It refers to the
difference between maximum possible complexity and the actual complexity and
in this case work means doing directed computation i.e. a computation with a
goal, which is called computational work unco2 .

In this paper, in order to study the anisotropy influence on HSC and provide
an informational view of the phase transition, we start to review an
anisotropic 5-dimensional background specified by an arbitrary dynamical
exponent, a nontrivial warp factor and two Maxwell fields that supports the
Van der Waals-like phase transition between small and large black holes, in
section II. Then we calculate HSC in this model for two cases of strip
entangling orientation, parallel and perpendicular to the isotropic
coordinate, in section III. Finally, we present our finding and numerical
results in section IV.

##  II Review on the background

Motivated by the experimental results about the dependence of total
multiplicity on energy, anisotropic holographic models related with the
Lifshitz-like background, have been studied ex1 ; ex2 . The metric we are
interested in is supported by a five-dimensional Einstein-Dilaton-two-Maxwell
action and introduced in metric  as follows

$ds^{2}=\frac{R^{2}b_{s}(z)}{z^{2}}\left(-g(z)dt^{2}+dx_{1}^{2}+z^{2-\frac{2}{\nu}}(dx_{2}^{2}+dx_{3}^{2})+\frac{dz^{2}}{g(z)}\right),$
(1)

where $z$ is radial direction and boundary coordinates where the field theory
lives on are $t,x_{1},x_{2}$ and $x_{3}$. It is anisotropic along the boundary
directions and we call $x_{1}$ and $\nu$, the isotropic coordinate and
anisotropic parameter, respectively. $R$ is the characteristic length scale of
the geometry that we set as equal to one in our numerical results.
$b(z)=e^{cz^{2}/2}$ is the warp factor. Since for $c<0$ the dilaton field
becomes real metric , in all our calculations we set $c=-1$. The blackening
function, $g(z)$, satisfies $g(z_{h})=0$ and $g(0)=1$ in which $z_{h}$ is the
position of the horizon, and it must hold the following equation metric :

$g^{\prime\prime}+g^{\prime}\left(\frac{3cz}{2}-\frac{1}{z}-\frac{2}{\nu
z}\right)-\frac{\mu^{2}c^{2}z^{2+\frac{2}{\nu}}e^{-cz^{2}}}{4\left(1-e^{-\frac{czh^{2}}{4}}\right)^{2}}=0,$
(2)

where $\mu$ is the chemical potential in the field theory, due to the
existence of the gauge field in the bulk metric . The Hawking temperature is
given as $T=|g^{\prime}(z_{h})|/4\pi$. We have plotted the temperature as a
function of horizon at different chemical potentials for isotropic (where
$\nu=1$) and anisotropic cases in Fig 1. In both cases, at $\mu=0$, the
temperature has a minimum $T_{\rm{min}}$ at $z_{h}=z_{\rm{min}}$. The branch,
corresponding to $0<z_{h}<z_{\rm{min}}$ (large black holes) is
thermodynamically stable and the branch in which temperature increases with
the horizon, $z_{h}>z_{\rm{min}}$, (small black holes) is unstable metric .
Below the minimum temperature there is no black hole solution and therefore,
there exists a Hawking-Page phase transition where, a thermal AdS solution
transits to the large black hole solution at temperature $T\geq T_{HP}$, and
$T_{HP}\mathchar 13350\relax T_{\rm{min}}$. For $0<\mu<\mu_{cr}$, the
temperature has a local minimum, $\bar{T}_{\rm{min}}$ and a local maximum,
$\bar{T}_{\rm{max}}$ at the values of the horizon, which we call
$\bar{z}_{\rm{min}}$ and $\bar{z}_{\rm{max}}$, respectively. The black holes
between $\bar{z}_{\rm{min}}$ and $\bar{z}_{\rm{max}}$ are thermodynamically
unstable and there is a black hole-black hole phase transition between small
black holes corresponding to $z_{h}>\bar{z}_{\rm{max}}$ and large black holes
corresponding to $0<z_{h}<\bar{z}_{\rm{min}}$ at $T_{BB}\mathchar
13350\relax\bar{T}_{\rm{min}}$ (metric, ). For $\mu>\mu_{cr}$, the temperature
decreases with the growth of the horizon, monotonously and there is no
transition anymore. In fact, a second order phase transition happens at the
critical point, i.e. $\mu=\mu_{cr}$. Note that as you can see in Fig 1, the
values of the critical point depend on the values of the anisotropic
parameter. To find more properties of this background, we refer the reader to
metric . To study the potential of quark pairs you can see metric2  and you
can find an investigation of holographic entanglement entropy in this
background in metric3 .

Figure 1: Temperature as a function of the horizon for isotropic case (left)
and anisotropic case, $\nu=2$ (right) for different values of $\mu$.

##  III Holographic subregion complexity

Inspired by the Hubney-Ryu-Takayanagi proposal, the CV proposal extends to
mixed states. HSC cojecture states that the complexity of the mixed state
defined on the boundary with a subregion like a rectangular region $A$ in our
case, is dual to the volume enclosed by the minimal hyper-surface
$\gamma_{A}$, appearing in the holographic entanglement entropy proposal
alishahiha , i.e.

${\cal{C}}_{A}=\frac{V_{\gamma_{A}}}{8\pi RG_{N}},$ (3)

where $R$ and $G_{5}$ are AdS radius and the five-dimensional Newton constant,
respectively. ${\cal{C}}_{A}$ is the HSC for the subregion $A$. The strip
entangling surface which is considered here, is defined with one dimension of
the length $l$ in $x_{1}$ in parallel to the isotropic coordinate, or in
$x_{2}$ or $x_{3}$ (since there is a rotational symmetry in these two
directions, there is no difference between them) in perpendicular to the
isotropic coordinate. In both cases the surface has a width
$L\rightarrow\infty$, see Fig 2. Using (1), we easily find the area of the
parallel and perpendicular minimal surfaces as follows

$\displaystyle
S^{||}=\frac{L^{2}}{4G_{5}}\int_{-l/2}^{l/2}\frac{R^{3}b(z)^{\frac{3}{2}}}{z^{3}}z^{2-\frac{2}{\nu}}\sqrt{1+\frac{z^{\prime}(x_{1})^{2}}{g(z)}}dx_{1},$
(4) $\displaystyle
S^{\bot}=\frac{L^{2}}{4G_{5}}\int_{-l/2}^{l/2}\frac{R^{3}b(z)^{\frac{3}{2}}}{z^{3}}z^{1-\frac{1}{\nu}}\sqrt{z^{2-\frac{2}{\nu}}+\frac{z^{\prime}(x_{2})^{2}}{g(z)}}dx_{2},$
(5)

where $z(x_{1})$ (or, equivalently, $x_{1}(z)$) in parallel case and
$z(x_{2})$ (or $x_{2}(z)$) in perpendicular case are the profile of the
minimal surface. By using the constant of motion, these profiles are obtained:

$\displaystyle
x_{1}(z)=2\int_{z}^{z_{*}}\frac{1}{\sqrt{g(z)}}\left(\frac{b(z)^{3}}{b(z_{*})^{3}}\frac{z_{*}^{2+\frac{4}{\nu}}}{z^{2+\frac{4}{\nu}}}-1\right)^{-\frac{1}{2}}dz,$
(6) $\displaystyle
x_{2}(z)=2\int_{z}^{z_{*}}\frac{1}{\sqrt{g(z)}}\left(\frac{b(z)^{3}}{b(z_{*})^{3}}\frac{z_{*}^{2+\frac{4}{\nu}}}{z^{\frac{6}{\nu}}}-z^{2-\frac{2}{\nu}}\right)^{-\frac{1}{2}}dz,$
(7)

where $z_{*}$ is the turning point of the minimal surface and lies at
$x_{1}=0$ in parallel case and $x_{2}=0$ in perpendicular case. In this static
background, by slicing the bulk with planes of constant $z$, the volume
enclosed by $\gamma_{A}$ is obtained by integrating the inside of the parallel
minimal surface as follows

Figure 2: Parallel minimal surface (green) and perpendicular minimal surface
(red).

$V(z_{*})^{||}=2L^{2}\int_{0}^{z_{*}}\frac{R^{4}b_{s}(z)^{2}}{z^{4}}z^{2-\frac{2}{\nu}}\sqrt{\frac{1}{g(z)}}x_{1}(z)dz.$
(8)

$V(z_{*})^{\bot}$ is obtained by replacing $x_{1}(z)$ with $x_{2}(z)$ in
equation (8). Putting anisotropic parameter $\nu=1$ in equation (8), leads to
$V_{iso}$. In this case, we have $x_{1}(z)=x_{2}(z)$ and
$V^{||}=V^{\bot}=V_{iso}$.

Volumes similar to the case of entanglement entropy are divergent. We
introduce a normalized volume using (3) for isotropic and anisotropic cases,
respectively,

$\displaystyle C_{iso}\equiv\frac{8\pi
RG_{5}}{L^{2}}(\mathcal{C}_{iso}-\mathcal{C}_{AdS})=\frac{V_{iso}-V_{AdS}}{L^{2}},$
(9) $\displaystyle C^{||or\perp}\equiv\frac{8\pi
RG_{5}}{L^{2}}(\mathcal{C}^{||or\perp}-\mathcal{C}_{0}^{||or\perp})=\frac{V^{||or\perp}-V_{0}^{||or\perp}}{L^{2}}.$
(10)

In equation (9), $\mathcal{C}_{iso}$ and $\mathcal{C}_{AdS}$ are the HSC for
$A$ in isotropic form of (1) and AdS background, respectively. The volumes are
defined for the same boundary region such that $V_{iso}$ reduces to $V_{AdS}$
by setting $c$ and $\mu$ equal to zero and $g(z)=1$. In equation (10),
$\mathcal{C}^{||or\perp}$ and $\mathcal{C}_{0}^{||or\perp}$ are the HSC for
$A$, parallel or perpendicular volumes in the anisotropic case, in the metric
(1) and anisotropic geometry with the same value of $\nu$ and $g(z)=1$ (zero
temperature solution), respectively.

##  IV Numerical results

In this section, we will present our findings from the numerical calculations
of normalized HSC, which is introduced in equations (9) and (10) in isotropic
and anisotropic cases.

###  IV.1 Isotropic background

In the isotropic background i.e. when we set the anisotropic parameter equal
to one, there is a rotational symmetry in three boundary spatial directions.
We would like to investigate whether the thermodynamically stable and unstable
solutions can be interpreted from a quantum information point of view. To
achieve this purpose, we study the behavior of $C_{iso}$ in terms of
temperature. In the left panel of Fig 3, $C_{iso}$ has been plotted as a
function of temperature for $\mu=0$. This figure shows that $C_{iso}$ is a
double-valued function, the large black holes branch and the small black holes
one, as expected. In comparison to the picture of temperature as a function of
horizon in the left panel of Fig 1, the dashed branch corresponds to the small
black holes (or thermodynamically unstable solutions) and the thick part
corresponds to the large black holes (thermodynamically stable solutions).
Interestingly, between two solutions with the same temperature, the solution
with smaller value of $C_{iso}$ belongs to the thermodynamically stable
branch. According to (9), the smaller value of $C_{iso}$ means that the amount
of information it takes to specify the considered state is less. Therefore, at
fixed temperature, the states described with less information, characterize
the stable branch. Note that $C_{iso}$, which is the difference between HSC
for the non-conformal background and conformal background is always negative,
indicating that non-conformal effects decrease HSC in this model.

Figure 3: left: $C_{iso}$ as a function of temperature for $\mu=0$ and
$l=0.5$. Right: The dependence of $C_{iso}$ on temperature at fixed $\mu=0.01$
and $l=0.5$.

It should be noted that there is a maximum value of HSC in the system we
investigate. In the left panel of Fig 3, this maximum value corresponds to the
value of $C_{iso}$ at $T_{min}\approx T_{HP}$, where a phase transition takes
place from the thermal AdS background to the large black hole solution. On the
other hand, considering the black hole as a computational machine, there is a
resource for performing a computation work, which is defined as the difference
between the maximum possible complexity and the actual complexity (unco2, ).
In the stable branch (thick line), with the increase of temperature, $C_{iso}$
deviates from this maximum value. This means that the actual complexity of the
stable branch tends to deviate from the maximum possible value of complexity
with the increase of temperature. Therefore, in the stable branch resource
increases with the increase of temperature and in other words we have a better
computational machine.

As shown in the right panel of Fig 3, at fixed $\mu=0.01$, between
$\bar{T}_{\rm{min}}$ and $\bar{T}_{\rm{max}}$, $C_{iso}$ is a multivalued
function of temperature with three branches. In the branch labeled with 1,
corresponding to the small black holes branch, $C_{iso}$ is constant with the
increase of temperature until $\bar{T}_{\rm{max}}$. The thermodynamically
unstable branch is displayed as a dashed line until $\bar{T}_{\rm{min}}$. The
second stable branch corresponding to the large black holes is labeled with 3.
Near $T_{BB}\mathchar 13350\relax\bar{T}_{\rm{min}}$, as temperature goes up,
phase transition happens from the small black holes solutions to the large
black holes solutions. It corresponds to a jump from the first stable branch
to the second one, $(1\rightarrow 3)$, corresponding to the smaller value of
$C_{iso}$. This means the states identified with less information are more
favorable thermodynamically. In short, in cases of $\mu=0$ and $\mu\neq 0$,
the state that is favorable thermodynamically, needs less information to be
specified from the quantum information perspective. Note that we started from
the small black hole branch and argued about what happens with rising
temperature.

Figure 4: $C_{iso}$ in terms of temperature at $l=0.5$ for three values of
$\mu$.

We want to consider the case that we do not have the assumption mentioned at
the end of the previous paragraph i.e we do not start from the small black
hole branch. In fact, we are interested in characterizing stable and unstable
solutions between the three branches. To reach this goal, we can use the
concept of resource. As is shown by the right panel of Fig 3, there are many
states with constant maximum value of HSC corresponding to the small black
hole branch. It seems that in the presence of $\mu$, as you can see in Fig 4,
if $\frac{\mu}{T}$ is small enough, there exist states whose HSC is maximum
and independent of the value of $\frac{\mu}{T}$. However, with the increase of
$\frac{\mu}{T}$, $C_{iso}$ depends on $\frac{\mu}{T}$. Due to the existence of
this maximum, we can apply the interpretation of a computational machine with
its resource for doing a computational work. Between these three branches the
smaller value of $C_{iso}$ is not enough to determine the stable and unstable
solutions. It is because, as shown in Fig 4, the branch labeled 2, which is an
unstable solution has less $C_{iso}$ than branch number 1. However, based on
definition of resource, in branch number 2 (the unstable branch), the actual
complexity is nearing the maximum possible value with the increase of
temperature. This means its resource as a computational machine, which is the
difference between actual complexity and maximum possible complexity,
decreases and thus we have a worse machine for computing with the increase of
temperature. Therefore, from the information theory point of view, the
solutions are (un) stable if with the increase of temperature, its resource
(decreases) increases. As it is summarized in table 1, zero and increasing
resource (branches number 1 and 3) corresponded to stable solutions, while
decreasing resource (branch number 2) characterizes the unstable solution.

In Fig 4, we have plotted $C_{iso}$ as a function of temperature for different
values of $\mu$. It shows that by nearing $\mu_{cr}$, the multivalued function
tends to be the unified one-value function. The dependence of $C_{iso}$ on
temperature at this unified branch corresponding to $\mu>\mu_{cr}$ is
presented.

Figure 5: Left: $C_{iso}$ in terms of $l$ for $\mu=0$ and four values of
temperatures near to $T_{HP}\approx 0.346$. Right: $C_{iso}$ as a function of
$l$ for $\mu=0.05$ and different values of the temperature below and above the
phase transition temperature, $T_{BB}\approx 0.344$.

The left panel of Fig 5 shows the behavior of $C_{iso}$ as a function of the
subregion length, $l$ for $\mu=0$. This figure has been plotted for the stable
branch. For small values of $l$, $C_{iso}$ is almost independent of
temperature. However, for larger values of $l$ the difference between values
of $C_{iso}$ grows with varying temperatures. In fact, with the increase of
$l$, the energy of the probe becomes comparable with $T_{HP}$ and $C_{iso}$
gives more information about the phase transition point. As this Fig shows, at
fixed $l$, by comparing non-conformal states, the closer to the phase
transition point, the more information needed to specify the state. Moreover,
nearing the phase transition temperature $T_{HP}$, this increase changes more
with varying temperatures.

In the right panel of Fig 5, we have plotted $C_{iso}$ in terms of $l$ at
$\mu=0.05$ for values of the temperature below (the small black hole branch)
and above (the large black hole branch) the phase transition temperature.
There is the same behavior with $\mu=0$ case for large black hole branches. As
you can see, by passing the phase transition point and jumping to the small
black hole branch, in a large enough $l$, there is a significant increase in
the value of $C_{iso}$ or equivalently, an increase in the needed information
to specify the non-conformal state. Another point is that in the small black
hole branch for large enough values of $l$, the probe reacts to nearing the
phase transition point by increasing the value of $C_{iso}$ and regardless of
how far or near the temperature is to the phase transition.

Table 1: A short summary of resource and stability $\rm{Branch}$ |  $\rm{Resource}$ |  $\rm{Stability}$
---|---|---
$1$ | $\rm{vanishes}$ | stable
$2$ | $\rm{decreasing}$ | unstable
$3$ | $\rm{increasing}$ | stable

###  IV.2 Anisotropic background

Here we consider the anisotropic background (1) for different values of $\nu$
in two cases, parallel and perpendicular to the isotropic coordinate as shown
in Fig 2. The numerical results have been presented as follows.

####  IV.2.1 Parallel case

Using (10), we have plotted $C^{||}$ as a function of temperature for
different values of $\mu$ in Fig 6. The left panel corresponding to the case
of $\mu=0$ shows the thermodynamically stable branch (thick branch) is
corresponding to the states with smaller value of $C^{||}$, i.e. ones that
need less information to be specified. Hence, like the isotropic case, the
states with less required information to be specified characterize the stable
branch. Again, there is a maximum value of HSC of the mixed state under study
at $T_{min}\approx T_{HP}$ and in the stable branch HSC deviates from this
maximum value with the increase of temperature. From the point of view of a
computational machine, it means that resource increases and we have a better
machine for computing with growing temperature. The left panel shows that with
the increase of $\nu$, the minimum of temperature and consequently $T_{HP}$
decreases.

Figure 6: Left: $C^{||}$ as a function of $T$ for three values of $\nu$ for
$l=0.5$ at $\mu=0$. Right: $C^{||}$ as a function of $T$ for three values of
$\mu$ at $l=0.5$ and $\nu=2$.

The right panel of Fig 6 for the case of $\mu\neq 0$, shows $C^{||}$ as a
function of $T$ is a multivalued function between corresponding
$\bar{T}_{min}$ and $\bar{T}_{max}$, which is shown in the right panel of Fig
1, for $\mu<\mu_{cr}$ in the anisotropic case. Passing $\mu_{cr}$, this
function tends to be a one-value function as expected. In this panel just like
the previous section, we can characterize the stable and unstable branches
from the information view point. It is noteworthy that as you can see from Fig
6, $C^{||}$ is negative. According to (10) it means the effect of temperature
on HSC is decreasing. In agreement with co15 , we can show that our previous
argument about this behavior based on an ensemble of microstates corresponding
to a given mixed macrostate is valid again. At zero temperature there exists
one microstate (unique configuration) but with increasing temperature, the
number of microstates corresponding to the mixed macrostate increases. Thus,
we specify the state with less information than we do at zero temperature.

Figure 7: In all three panels $l=0.5$ and in the isotropic case $\nu=1$ and in
the anisotropic case $\nu=1.01$. $\Delta C^{||}$ as a function of $T$ at
$\mu=0$ for the stable branch, that is large black holes ones (left). $\Delta
C^{||}$ as a function of $T$ at $\mu=0.05$ for the stable branches i.e. large
black holes (middle) and small black holes ones (right).

In order to study the effect of anisotropic parameter on HSC in our system, we
have plotted $\Delta C^{||}\equiv\frac{8\pi
RG_{5}}{L^{2}}(\mathcal{C}^{||}-\mathcal{C}_{iso})$ as a function of
temperature for stable solutions in Fig 7. The left panel is $\Delta C^{||}$
for stable branch at $\mu=0$, and the middle and right panel are related to
the case of $\mu\neq 0$ for stable branches, large black holes (middle) and
small black holes (right). It shows $\Delta C^{||}$ is negative which means
anisotropy causes decrease of HSC. We want to justify this decrease by saying
that the system we study consists of a state and its environment which forms a
closed system. In this case, the external source that injects anisotropy into
the boundary field theory is the environment. In the initial situation the
source of anisotropy is zero and therefore the information for specifying it
is zero too. Thus, the needed information for specifying the whole closed
system is only related to the state. However, during the assumed evolution, in
the final situation, the external source value becomes non-zero and we need
information to specify it. According to the principle of information
conservation which states information is never lost from a closed isolated
system information , we claim with increase of complexity of the external
source in the final situation, the complexity of the state must decrease, as
you can see from Fig 7. We have proposed the same justification for an
external heat source in co15 .

Figure 8: Left: $C^{||}$ as a function of $l$ at $\nu=2$ and $\mu=0$ for
temperatures above and near $T_{HP}\approx 0.290$. Right: $C^{||}$ as a
function of $l$ at $\nu=2$ for $\mu=0.05$ and different values of temperature
below and above the phase transition temperature, $T_{BB}\approx 0.288$.

Fig 8 has presented the dependence of $C^{||}$ on subregion length at
different values of temperature. The left panel shows $C^{||}$ as a function
of $l$ at $\mu=0$ and has been plotted for large black holes i.e. the stable
branch. $C^{||}$ is almost independent of temperature for small values of $l$,
however with the increase of $l$ the difference is revealed, which compared
with the isotropic case, the left panel of Fig 5, this takes place at smaller
$l$. The right panel of Fig 8 shows $C^{||}$ as a function of $l$ at
$\mu=0.05$ for temperatures below (the small black hole branch) and above (the
large black hole branch) the phase transition temperature, $T_{BB}$. In this
case, $C^{||}$ undergoes a jump when the temperature crosses the $T_{BB}$.

####  IV.2.2 Perpendicular case

Using (10) for the perpendicular case, we have plotted $C^{\bot}$ in terms of
temperature in Fig 9. The left panel shows the case of $\mu=0$ in which
$C^{\bot}$ as a function of temperature is a double-valued function for three
values of $\nu$ as expected. The same argument of the previous subsections
applies here, for specifying the stable and unstable branch. The right panel
presents $C^{\bot}$ in terms of temperature as a multivalued function for
$\mu<\mu_{cr}$ and a one-value function by passing $\mu_{cr}$. To characterize
the stable branches and unstable branches we can use the same arguments as
before.

Figure 9: Left: $C^{\bot}$ as a function of $T$ for three values of $\nu$ for
$l=0.5$ at $\mu=0$. Right: $C^{\bot}$ in terms of $T$ for different values of
$\mu$ at $l=0.5$ and $\nu=2$.

To study the effect of the anisotropic parameter on HSC in this case, we have
plotted $\Delta C^{\bot}\equiv\frac{8\pi
RG_{5}}{L^{2}}(\mathcal{C}^{\bot}-\mathcal{C}_{iso})$ as a function of
temperature for stable solutions in Fig 10. The left panel indicates $\Delta
C^{\bot}$ at $\mu=0$ case for stable branch (large black holes), and the
middle and right panel shows the case of $\mu\neq 0$ for stable branches,
large black holes (middle) and small black holes (right). Like the parallel
case, the anisotropy effect on HSC of the state under probe is a decreasing
effect. We justify this reduction by considering the whole closed system
consists of a state and its environment, just like what was described in
detail in the previous section for the parallel case. Another thing that is
clear from this Fig compared to Fig 7 is that $C^{||}>C^{\bot}$. We have no
argument to justify it at this point, but we guess it might be related to the
screening effect. In the sense that the effects of anisotropy causes a
stronger screening in the perpendicular direction and thus we need less
information to specify the state in this case, in agreement with co14 .

Figure 10: In all three panels $l=0.5$ and in the isotropic case $\nu=1$ and
in the anisotropic case $\nu=1.01$. $\Delta C^{\bot}$ in terms of $T$ at
$\mu=0$ for the stable branch, that is large black holes ones (left). $\Delta
C^{\bot}$ as a function of $T$ at $\mu=0.05$ for the stable branches i.e.
large black holes (middle) and small black holes ones (right).

Fig 11 indicates $C^{\bot}$ as a function of $l$ for different values of
temperature. The left panel is the case of $\mu=0$ and has been plotted for
stable branch (large black holes) which approach the phase transition
temperature, $T_{HP}$ from above. The difference between values of $C^{\bot}$
appears with the increase of $l$, which compared with the isotropic case, this
takes place at smaller $l$, like the parallel case. The right panel shows
$C^{\bot}$ as a function of $l$ for $\mu=0.05$. This panel has been plotted
for temperatures below (the small black holes) and above (the large black
holes) the phase transition temperature, $T_{BB}$. In this case, like parallel
case, $C^{\bot}$ has a jump when the temperature passes the $T_{BB}$.

Figure 11: Left: $C^{\bot}$ as a function of $l$ at $\nu=2$ and $\mu=0$ for
temperatures above and near $T_{HP}\approx 0.290$. Right: $C^{\bot}$ as a
function of $l$ at $\nu=2$ for $\mu=0.05$ and different values of temperature
below and above the phase transition temperature, $T_{BB}\approx 0.288$.

## References

  * (1) J. Casalderrey-Solana, H. Liu, D. Mateos, K. Rajagopal and U. A. Wiedemann, “Gauge/String Duality, Hot QCD and Heavy Ion Collisions,” Cambridge University Press, 2014. [arXiv:1101.0618 [hep-th]].
  * (2) R. Baier, A. H. Mueller, D. Schiff and D. T. Son, “Bottom up thermalization in heavy ion collisions,” Phys. Lett. B 502, 51 (2001) [arXiv:0009237 [hep-th]].
  * (3) M. Strickland, “Thermalization and isotropization in heavy-ion collisions,” Pramana 84, 671, (2015) [arXiv:1312.2285 [hep-th]].
  * (4) D. Mateos and D. Trancanelli, “The anisotropic N=4 super Yang-Mills plasma and its instabilities,” Phys.Rev.Lett 107, 101601, (2011) [arXiv:1105.3472 [hep-th]].
  * (5) D. Giataganas and U. Gursoy and J. F. Pedraza, “Strongly-coupled anisotropic gauge theories and holography,” Phys.Rev.Lett 121, 121601, (2018) [arXiv:1708.05691 [hep-th]].
  * (6) I.Y. Aref’eva, A.A. Golubtsova, “Shock waves in Lifshitz-like spacetimes,” JHEP 1504, 011 (2015) [arXiv:1410.4595 [hep-th]].
  * (7) Micheal A. Nielsen and Isaac L. Chuang, “Quantum Computation and Quantum Infromation,” Cambridge university press, (2010), 702 p.
  * (8) R. Jefferson and R. C. Myers “Circuit complexity in quantum field theory,” JHEP 10 (2017) 107, [arXiv:1707.08570[quant-ph]].
  * (9) D. Stanford and L. Susskind, “Complexity and Shock Wave Geometries,” Phys. Rev. D 90, no.12, 126007 (2014) [arXiv:1406.2678[hep-th]].
  * (10) A. R. Brown, D. A. Roberts, L. Susskind, B. Swingle and Y. Zhao, “Complexity, action, and black holes,” Phys. Rev. D 93, no.8, 086006 (2016) [arXiv:1512.04993[hep-th]].
  * (11) Mohsen Alishahiha, “Holographic complexity,” Phys. Rev. D. 92, 126009 (2015) [arXiv:1509.06614[hep-th]].
  * (12) Omer Ben-Ami, Dean Carmi, “On Volumes of Subregions in Holography and Complexity,” JHEP 1611, 129 (2016) [arXiv:1609.02514[hep-th]].
  * (13) S. J. Zhang, “Complexity and phase transitions in a holographic QCD model,” Nucl. Phys. B 929, 243 (2018) [arXiv:1712.07583[hep-th]].
  * (14) S. J. Zhang, “Subregion complexity in holographic thermalization with dS boundary,” Eur.Phys.J.C 79 (2019) 8,715 [arXiv:1905.10605[hep-th]].
  * (15) Pratim Roy, Tapobrata Sarkar, “On subregion holographic complexity and renormalization group flows,” Phys.Rev. D 97, 086018 (2018) [arXiv:1708.05313[hep-th]].
  * (16) R. Fareghbal and P. Karimi, “Complexity growth in flat spacetimes,” Phys. Rev. D 98, no. 4, 046003 (2018) [arXiv:1806.07273[hep-th]].
  * (17) M.Alishahiha, A.Faraji Astaneh, M.R.Mohammadi Mozaffar and A.Mollabashi, “Complexity Growth with Lifshitz Scaling and Hyperscaling Violation,” JHEP 1807, 042 (2018) [arXiv:1802.06740 [hep-th]].
  * (18) M. Alishahiha, K. Babaei Velni and M. R. Mohammadi Mozaffar, “Subregion Action and Complexity,” Phys.Rev.D 99 (2019) 12, 126016, [arXiv:1809.06031 [hep-th]].
  * (19) M. asadi, “On volume subregion complexity in non-conformal theories,” Eur.Phys.J.C 80 (2020) 7, 681 [arXiv:2004.11306[hep-th]].
  * (20) Mahsa Lezgi, Mohammad Ali-Akbari and Mohammad Asadi, “Non-Conformality, Subregion Complexity and Meson Binding,” Phys.Rev.D 104, 026001 (2021) [arXiv:2011.11625[hep-th]].
  * (21) Mahsa Lezgi and Mohammad Ali-Akbari “Note on holographic subregion complexity and QCD phase transition, ” Phys. Rev. D 101, 026022 (2020) [arXiv:1908.01303[hep-th]].
  * (22) Mahsa lezgi and Mohammad Ali-Akbari, “Complexity and uncomplexity during energy injection,” Phys. Rev. D 103 (2021) 12, 126024 [arXiv:2103.05023[hep-th]].
  * (23) Leonard Susskind, “Three Lectures on Complexity and Black Holes,” [arXiv:1810.11563[hep-th]].
  * (24) Adam R. Brown and Leonard Susskind, “Second law of quantum complexity,” Phys. Rev. D 97 (2018) 8, 086015 [arXiv:1701.01107[hep-th]].
  * (25) Irina Aref’eva, Anastasia A. Golubtsova, Eric Gourgoulhon “Analytic black branes in Lifshitz-like backgrounds and thermalization,” JHEP 09, 142 (2016) [arXiv:1601.06046[hep-th]].
  * (26) Irina Aref’eva, “Holography for Heavy Ions Collisions at LHC and NICA,” EPJ Web Conf 164, 01014 (2017) [arXiv:1612.08928[hep-th]].
  * (27) Irina Aref’eva, Kristina Rannu, “Holographic Anisotropic Background with Confinement-Deconfinement Phase Transition,” JHEP 05, 206 (2018) [arXiv:1802.05652[hep-th]].
  * (28) Irina Aref’eva, Kristina Rannu and Pavel Slepov, “Orientation Dependence of Confinement-Deconfinement Phase Transition in Anisotropic Media,” Phys.Lett.B 792, 470-475 (2019) [arXiv:1808.05596[hep-th]].
  * (29) Irina Aref’eva, Alexander Patrushev and Pavel Slepov, “Holographic entanglement entropy in anisotropic background with confinement-deconfinement phase transition,” JHEP 07, 043 (2020) [arXiv:2003.05847[hep-th]].
  * (30) L. Susskind and J. Lindesay, “An introduction to black holes, information and the string theory revolution: The holographic universe,” Hackensack, USA: World Scientific (2005) 183 p.


# The arions generation by magnetodipole waves
of pulsars and magnetars in a constant magnetic field

V.I.Denisov1, G.A.Dantsev1, V.I.Priclonsky1,
I.P.Denisova2,3,a, O.N.Gavrish2,3
1Department of Physics, Moscow State University,
119991, Moscow, Russia,
2 Bauman Moscow State Technical University,
105005, Moscow, Russia
3 The State University of Management,
109542, Moscow, Russia
a e-mail<EMAIL_ADDRESS>(corresponding author)

The influence of the gravitational fields of pulsars and magnetars on the
arion emission during the propagation of magnetodipole waves in a constant
magnetic field has been evaluated.

The solution of the equation was obtained and the flux of arions emitted by
magnetodipole waves during their propagation in a constant magnetic field was
found. It is shown that the amplitude of the born arion wave at a distance
from the source of magnetodipole radiation of a pulsar or magnetar
$(r\to\infty)$ in the considered case tends to a constant value. The intensity
of the arion emission in the solid angle element and the amount of arion
energy $\overline{I}$, emitted in all directions per unit time grow
quadratically with increasing distance, traveled by the magnetodipole
radiation of a pulsar or magnetar in a constant magnetic field.

Such growth of the energy of the born arion wave is due to the fact that in
the considered problem constant magnetic field is defined in the whole space.
In reality, the galactic and intergalactic magnetic fields can be represented
in this form only in regions of space of finite dimensions, outside of which
the force lines of their induction vector are curved. Therefore, it is
possible to apply these results only in a region of space for which $r\leq
L_{coh}<\infty$, where $L_{coh}$ is the coherence length, the distance at
which the force lines of the induction vector can be considered as straight
lines. An estimate for the value of the coupling constant of photons with
arions is obtained.

## 1\. Introduction

In the scientific literature of recent years, the processes of photoproduction
of various axion-like particles beyond the Standart Model: arions [1,2] axions
[3-6], and dilatons [7-10] are actively discussed. These processes are
currently regarded as the most realistic processes, with the help of which it
is supposed [11-13] to carry out registration of axion-like particles in
laboratory and astrophysical conditions.

The arion is a strictly massless pseudoscalar Goldstone particle $a,$ which
was introduced in 1982, in the papers [14-16] of Prof. A. A. Anselm and his
co-authors.

The density of the Lagrangian function for the arion field, which is
interacting with the electromagnetic field, is usually written in the
canonical form:

$L={\sqrt{-g}\over 2}g^{nm}{\partial a\over\partial x^{n}}{\partial
a\over\partial x^{m}}-{\sqrt{-g}\over
16\pi}F_{nm}F^{nm}-{g_{a\gamma}\sqrt{-g}\over 4}F_{nm}\tilde{F}^{nm}a,$ $None$

where $a$ is the pseudoscalar field of the arion, $g_{a\gamma}$– coupling
constant of the arion with the electromagnetic field, $F_{nm}$ –
electromagnetic field tensor, $g$ is the determinant of the metric tensor,
$\tilde{F}^{nm}=E^{nmik}F_{ik}/2$ and $E^{nmik}=e^{nmik}/\sqrt{-g}$ – the
axial absolutely antisymmetric Levi-Civita tensor, and $e^{nmik}$ is the axial
absolutely antisymmetric Levi-Civita symbol, and $e^{0123}=+1.$

When studying the processes of arion generation under astrophysical conditions
the influence of the gravitational field, in general, cannot be neglected.
Therefore, first of all, let us estimate the magnitude of this influence and
the size of the region of space in which this influence can be significant.
Since the distribution of matter in neutron stars is close to spherically
symmetric, then we will use the Schwarzschild solution as the metric tensor of
the pseudo-Riemannian space. In the paper [18] it is shown, that the
Schwarzschild solution can be the external solution for a non-spherically
symmetric distribution of matter.

The most convenient coordinates for writing down this solution are isotropic
coordinates [17]. The nonzero components of the metric tensor in these
coordinates have the form:

$g_{00}={(4r-r_{g})^{2}\over(4r+r_{g})^{2}},\ \ \ \ \ \
g_{xx}=g_{yy}=g_{zz}=-(1+{r_{g}\over 4r})^{4},$ $None$

where $r_{g}$ is the Schwarzschild radius and the notation is used for
convenience of reference: $r=\sqrt{x^{2}+y^{2}+z^{2}}.$

Let us estimate the value of the ratio $r_{g}/r$, included in expressions (2),
on the surface of a neutron star. The radius of a neutron star in recent times
[19] is taken to be $R_{s}=10$ km, and the mass $M$ in the interval from 0.1
to 1.0 solar masses.

Thus, in our problem, the ratio $r_{g}/R_{s}\sim 0.01$. Since at $r>R_{s}$ the
ratio $r_{g}/r$ takes smaller values, than the ratio $r_{g}/R_{s}$, the
influence of the gravitational field is small and limited to a small
neighborhood $r\leq 100R_{s}$ of the neutron star. Therefore, as a first
approximation for the small parameter $r_{g}/r$ in our problem as a metric
tensor we will use the metric tensor of the Minkowski space: $g_{00}=1,\
g_{11}=g_{22}=g_{33}=-1.$

The equations of the arion and electromagnetic fields derived from the
Lagrangian density (1), in the Minkowski space have the form:

$\vbox{\hrule\hbox to5.0pt{\vrule width=0.5pt,height=5.0pt\hfill\vrule
width=0.5pt,height=5.0pt} \hrule}\
a=-{g_{a\gamma}}F_{mp}\tilde{F}^{mp}=-{g_{a\gamma}}({\bf B\ E}),$ $None$
${\partial F^{nm}\over\partial x^{m}}=-4\pi\tilde{F}^{nm}{\partial
a\over\partial x^{m}},$

where ${\bf B}$ is the magnetic field induction vector, ${\bf E}$ is the
electric field strength vector.

According to the first equation (3), the source of arions are electromagnetic
fields and waves, in which the first invariant of the electromagnetic field
tensor is different from zero. In the nature there exist such configurations
of electromagnetic fields and waves at which this invariant is different from
zero in large, by earthly standards, volumes of space. These are, for example,
magnetodipole radiation of pulsars and magnetars, propagating in the constant
galactic or intergalactic magnetic field. And although the induction of
galactic and intergalactic magnetic fields is relatively small at $B\sim
10^{-6}$ Gs, the volumes occupied by these fields are significant. Therefore,
it is of undoubted interest to study this process. Let us consider it in
detail.

## 2\. Calculation of arion emission arising from the propagation of
magnetodipole waves of a pulsar or a magnetar in a constant magnetic field

Suppose that at a point with radius vector ${\bf r}={\bf
r}_{0}=\\{x_{0},y_{0},z_{0}\\}$ is a pulsar or magnetar with a magnetic dipole
moment $\bf m$, which rotates at a frequency of $\omega$ around an axis that
makes an angle of $\psi$ with the vector $\bf m$. Then this source emits an
electromagnetic wave of frequency $\omega$, whose components [8] have the
form:

${\bf B}({\bf R},\tau)={3({\bf m}(\tau)\cdot{\bf R}){\bf R}-R^{2}{\bf
m}(\tau)\over R^{5}}-{{\dot{\bf m}}(\tau)\over cR^{2}}+$ $None$ $+{3({\dot{\bf
m}}(\tau)\cdot{\bf R}){\bf R}\over cR^{4}}+{(\ddot{\bf m}(\tau)\cdot{\bf
R}){\bf R}-R^{2}\ddot{\bf m}(\tau)\over c^{2}R^{3}},$ ${\bf E}({\bf
R},\tau)={({\bf R}\times\dot{\bf m}(\tau))\over cR^{3}}+{({\bf
R}\times\ddot{\bf m}(\tau))\over c^{2}R^{2}},$

where ${\bf R}={\bf r}-{\bf r}_{0},$ the dot above the vector means the
derivative on retarded time $\tau=t-R/c$, and the pulsar or the magnetar
magnetic dipole moment in the following task has the components:

${\bf m}(\tau)=|{\bf m}|\\{\cos(\omega\tau)\sin\psi,\
\sin(\omega\tau)\sin\psi,\ \cos\psi\\}.$ $None$

In the coordinate system, the origin of which is placed in the point ${\bf
r}={\bf r}_{0}$, the directional diagram of electromagnetic radiation (4) has
the form

$\overline{\left({dI\over d\Omega}\right)}_{EMW}={[{\bf r}\ddot{\bf
m}]^{2}\over 4\pi c^{3}r^{2}}.$

Integrating this expression over the angles $\theta$ and $\varphi$, we obtain
the total intensity of radiation of a pulsar or a magnetar:

$\overline{I}_{EMW}={2\ddot{\bf m}^{2}\over 3c^{3}}={2ck^{4}{\bf
m}^{2}\sin^{2}\psi\over 3}.$ $None$

As shown in [2,8], the electromagnetic wave (4) can emit arions and dilatons
with frequencies $\omega$ and $2\omega$.

Let there is a constant magnetic field in the considered region of space

${\bf B}_{0}=\\{B_{0x},B_{0y},B_{0z}\\}.$ $None$

Substituting expressions (4) and (7) into equation (3) and considering that

${\bf E}({\bf r},\tau)={({\bf R}\times\dot{\bf m}(\tau))\over cR^{3}}+{({\bf
R}\times\ddot{\bf m}(\tau))\over c^{2}R^{2}}=rot_{{\bf r}_{0}}{\dot{\bf
m}(\tau)\over cR},$

where $rot_{{\bf r}_{0}}$ is the rotor taking operator on coordinates of the
vector ${\bf r}_{0}=\\{x_{0},y_{0},z_{0}\\}$, let’s put it in the form

$\vbox{\hrule\hbox to5.0pt{\vrule width=0.5pt,height=5.0pt\hfill\vrule
width=0.5pt,height=5.0pt} \hrule}\ a=-{g_{a\gamma}}({\bf B}_{0}\ rot_{{\bf
r}_{0}}{\dot{\bf m}(\tau)\over cR}).$ $None$

Now let’s take into account that

$({\bf B}_{0}\ rot_{{\bf r}_{0}}{\dot{\bf m}(\tau)\over cR})=div_{{\bf
r}_{0}}\left[{\dot{\bf m}(\tau)\over cR}\ {\bf B}_{0}\right].$

Introducing an auxiliary vector ${\bf F}$, we express the pseudoscalar field
$a$ through it:

$a=div_{{\bf r}_{0}}\left[{\bf F\ B}_{0}\right].$

Then, to fulfill equation (8), it is necessary to require that this vector
satisfies the equation:

$\vbox{\hrule\hbox to5.0pt{\vrule width=0.5pt,height=5.0pt\hfill\vrule
width=0.5pt,height=5.0pt} \hrule}\ {\bf F}=-{g_{a\gamma}}{\dot{\bf
m}(\tau)\over cR}.$ $None$

Let’s substitute expression (5) into the right part of this equation and for
convenience of its solution we will write equation (9) in complex form:

$\vbox{\hrule\hbox to5.0pt{\vrule width=0.5pt,height=5.0pt\hfill\vrule
width=0.5pt,height=5.0pt} \hrule}\ {\bf F}=-{{g_{a\gamma}}|{\bf
m}|\omega\sin\psi\over c}\Big{(}i{\bf e}_{x}+{\bf e}_{y}\Big{)}{\exp{i(\omega
t-kR)}\over R}.$ $None$

After solving this equation, let’s take only the real part.

We will look for the solution of equation (10) in the form:

${\bf F}=-{{g_{a\gamma}}|{\bf m}|\omega\sin\psi\over c}\Big{(}i{\bf
e}_{x}+{\bf e}_{y}\Big{)}Q,$

where the function $Q$ must satisfy the equation:

$\vbox{\hrule\hbox to5.0pt{\vrule width=0.5pt,height=5.0pt\hfill\vrule
width=0.5pt,height=5.0pt} \hrule}\ Q={\exp{i(\omega t-kR)}\over R}.$

The retarded solution of this equation is of the form:

$Q={i\over 2k}\exp{i(\omega t-kR)}.$

As a result, for the real part of axion radiation arising at propagation of
the magnetodipole radiation (4) of a pulsar or magnetar in a constant magnetic
field (7), we have:

$a={g_{a\gamma}|{\bf m}|k\sin\psi\over
2R}\Big{\\{}B_{0z}\Big{[}(x-x_{0})\cos[\omega t-kR]+(y-y_{0})\sin[\omega
t-kR]\Big{]}-$ $None$ $-(z-z_{0})\Big{[}B_{0y}\sin[\omega
t-kR]+B_{0x}\cos[\omega t-kR]\Big{]}\Big{\\}}.$

Let us pass to the Cartesian frame of reference, the origin of which is placed
in the center of the pulsar or magnetar. Then $x_{0}=y_{0}=z_{0}=0$ and
expression (11) will take the form:

$a={g_{a\gamma}|{\bf m}|k\sin\psi\over 2r}\Big{\\{}B_{0z}\Big{[}x\cos[\omega
t-kr]+y\sin[\omega t-kr]\Big{]}-$ $None$ $-z\Big{[}B_{0y}\sin[\omega
t-kr]+B_{0x}\cos[\omega t-kr]\Big{]}\Big{\\}}.$

Thus, the arion emission arising from the propagation of the magnetodipole
radiation of a pulsar or magnetar radiation in a constant magnetic field
occurs at the same frequency, as the frequency of the electric field that
generates it.

## 3\. Angular distribution of arion radiation

Angular distribution of the arion emission, arising from the propagation of
the electromagnetic wave (4) of a pulsar or magnetar through a permanent
magnetic field (7), following the paper [8], we write it in the form:

${dI\over d\Omega}=r\big{(}{\bf r\ W}\big{)},$ $None$

where $\bf W$ is the energy flux density vector associated with the components
of the energy-momentum tensor $T^{nk}$ by the relation
$W^{\beta}=cT^{0\beta}$. Using the expression for the energy-momentum tensor
of the free arion field

$T^{nk}=g^{np}g^{km}\big{\\{}{\partial a\over\partial x^{p}}{\partial
a\over\partial x^{m}}-{1\over 2}g^{nk}{\partial a\over\partial x^{p}}{\partial
a\over\partial x^{m}}g^{pm}\big{\\}},$

from expression (13) we obtain:

${dI\over d\Omega}=-cr^{2}{\partial a\over\partial r}{\partial a\over\partial
x^{0}}.$

Substituting into this relation the expression (12) for the arion field and
leaving asymptotically the main part, after time averaging, we reduce it to
the form:

$\overline{{dI\over d\Omega}}={g^{2}_{a\gamma}c|{\bf
m}|^{2}k^{4}r^{2}\sin^{2}\psi\over
8}\Big{\\{}(B^{2}_{0x}+B^{2}_{0y})\cos^{2}\theta+B^{2}_{0z}\sin^{2}\theta-$
$None$
$-2B_{0z}\sin\theta\cos\theta[B_{0x}\cos\varphi+B_{0y}\sin\varphi]\Big{\\}}.$

Let us now find the amount of arion energy $\overline{I}$ emitted in all
directions per unit time:

${\overline{I}}_{AR}=\int\limits_{0}^{\pi}\sin\theta
d\theta\int\limits_{0}^{2\pi}d\varphi{\overline{dI}\over d\Omega}={\pi
g^{2}_{a\gamma}c|{\bf m}|^{2}k^{4}r^{2}\sin^{2}\psi\over
12}\Big{\\{}B^{2}_{0x}+B^{2}_{0y}+2B^{2}_{0z}\Big{\\}}.$ $None$

The different dependence of this result on the components of the constant
magnetic field (7) is due to the fact that the directivity diagram of the
magnetodipole radiation of a pulsar or a magnetar is not spherically
symmetric.

## 4\. Conclusion

It follows from expression (12) that the amplitude of the born arion wave at a
distance from the magnetodipole source of the pulsar or magnetar
$(r\to\infty)$ in the considered case tends to a constant value. The intensity
of the arion emission in the solid angle element and the amount of arion
energy $\overline{I}$, emitted in all directions per unit time grow
quadratically with increasing distance, traveled by the magnetodipole
radiation of a pulsar or magnetar in a constant magnetic field (7).

Such growth of the energy of the born arion wave is due to the fact that in
our problem the constant magnetic field (7) is set in the whole space. In
reality, the galactic and intergalactic magnetic fields can be represented in
the form (7) only in regions of finite dimensions, outside of which the force
lines of their induction vectors are curved. Therefore, one can apply the
results (12), (14) and (15) only in the region of space, for which $r\leq
L_{coh}<\infty$, where $L_{coh}$ is the coherence length – the distance at
which the lines of force of the induction vector can be considered as straight
lines.

Let us define the conversion factor of electromagnetic radiation energy into
arion energy $\beta$ as the ratio of the intensity of radiation of arions
${\overline{I}}_{ar}$ to the intensity of the electromagnetic field of a
pulsar or magnetar ${\overline{I}}_{EMW}$:

$\beta={{\overline{I}}_{AR}\over{\overline{I}}_{EMW}}={\pi
L^{2}_{cog}g^{2}_{a\gamma}\over
8}\Big{\\{}B^{2}_{0x}+B^{2}_{0y}+2B^{2}_{0z}\Big{\\}}.$

This coefficient reflects the properties of the converter, i.e., the
electromagnetic field. It should be noted, that according to the second
equation of the system (3), along with the conversion of the energy of
electromagnetic waves into the energy of arions, the opposite process takes
place. If $\beta<<1,$ then the reverse process can be ignored; if the
coefficient of $\beta$ is close to unity, it is necessary to consider both
processes together.

At present the value of the photon-arion coupling constant $g_{a\gamma}$ is
unknown. Let’s roughly estimate the value of this constant using the example
of the generation of arions by the electric field of pulsar radiation as it
propagates through the magnetic field of our Galaxy.

According to modern data [20,21], the radius of the Galaxy is 16 kiloparsecs.
It is accepted to distinguish the large-scale component of the magnetic field
of the Galaxy magnetic field of the Galaxy (the scale of homogeneity of the
order of hundreds and thousands of parsecs) and fluctuation component with a
wide range of scales (from fractions of parsecs to hundreds of parsecs). The
induction of the large-scale magnetic field of the Galaxy is estimated to be
2-3 $\mu$ Gs.

Based on these data, it is reasonable to put $L_{coh}\sim 10^{3}$ ps $=3\cdot
10^{21}$ cm, $B_{0}\sim 10^{-6}$ Gs.

Then from the condition $\beta<<1$ we obtain: $g^{2}_{a\gamma}<8\cdot
10^{-31}$ ${cm\over erg}.$ Equation (15) was derived by us into the Gaussian
system of units of measurement. Let us rewrite this formula in the natural
system of units. If we take into account that 1 erg= 624 GeV, and 1 cm =
$0.5\cdot 10^{14}$ GeV-1, then we get: $g_{a\gamma}<0.9\cdot 10^{-10}$ GeV-1.
This estimate coincides with the estimates of the coupling constant
$g_{a\gamma}$ obtained in [2,22-25].

## Acknowledgements

This study was conducted within the scientific program of the National Center
for Physics and Mathematics, section $\\#5$ $<<$Particle Physics and
Cosmology$>>.$ Stage 2023-2025.

## References

  * [1] A.A. Anselm, Experimental test for arion-photon oscillations in a homogeneous constant magnetic field, Phys. Rev. D 37, 2001 (1988)

  * [2] V.I. Denisov, B.D. Garmaev, I.P. Denisova, Radiation of arions by electromagnetic field of rotating magnetic dipole, Phys.Rev. D 104, 055018 (2021)

  * [3] N. Okada, D. Raut, Q. Shafi, Axions, WIMPs, proton decay and observable r in SO(10), Eur. Phys. J. C 83, 273 (2023)

  * [4] G. Raffelt, Astrophysical methods to constrain axions and other novel particle phenomena, Phys. Rep. 198, 1 (1990)

  * [5] R. Casadio, R. da Rocha, Axion stars in MGD background, Eur. Phys. J. C 83, 537 (2023)

  * [6] Y. Giraldo, R. Martinez, E. Rojas, J. Salazar, A minimal axion model for mass matrices with five texture-zeros, Eur. Phys. J. C 83, 638 (2023)

  * [7] I.P. Denisova, Dilatons generation while plane electromagnetic wave propagates in Coulomb field, Gravit. Cosmol. 27, 392 (2021)

  * [8] V.I. Denisov, I.P. Denisova, E.T. Einiev, The investigation of lowfrequency dilaton generation, Eur. Phys. J. C 82, 311 (2022)

  * [9] A. Ahmed, S. Najjari, Ultraviolet freeze-in dark matter through the dilaton portal, Phys. Rev. D 107, 055020 (2023)

  * [10] M. Astashenkov, Dilaton photoproduction in a magnetic dipole field of pulsars and magnetars, Eur. Phys. J. C 83 643 (2023)

  * [11] D. Banerjee et al. (NA64 Collaboration), Search for axionlike and scalar particles with the NA64 experiment, Phys. Rev. Lett. 125, 081801 (2020)

  * [12] Pierluca Carenzaa , Pedro De la Torre Luqueb, Detecting neutrino-boosted axion dark matter in the MeV gap, Eur. Phys. J. C 83 110 (2023)

  * [13] M. Hajebrahimi, H. Manshouri, M. Sharifian, M. Zarei, Axion-like dark matter detection using Stern-Gerlach Interferometer, Eur. Phys. J. C 83 11 (2023)

  * [14] A.A. Anselm, N. G. Uraltsev, A second massless axion?, Phys. Lett. 114B, 39 (1982)

  * [15] A.A. Anselm, Possible new long-range interaction and methods for detecting it, JETP Lett. 36, 55 (1982)

  * [16] A.A. Anselm, N.G. Uraltsev, Long range ”arion” field in the radio frequency band, Phys. Lett. 116B, 161 (1982)

  * [17] L. D. Landau, E. M. Lifshitz, The Classical Theory of Fields (Butterworth -Heinemann, Washington, DC, 1975)

  * [18] J.L. Hernandez-Pastora, L. Herrera, Non-spherical sources of Schwarzschild space-time, Eur. Phys. J. C 83 375 (2023)

  * [19] C.D. Capano, I. Tews, S.M. Brown, B. Margalit, S. De, S. Kumar, D.A. Brown, B. Krishnan, S. Reddy, Stringent constraints on neutron-star radii from multimessenger observations and nuclear theory, Nat. Astron 4, 625 (2020)

  * [20] V.G.Surdin. Distance to the galactic centre, Astron. Astrophys. Trans. 18, 2 367-372 (1999)

  * [21] F.Combes, J.Lequeux. The Milky Way: Structure, Dynamics, Formation and Evolution, Les Ulis: EDP Sciences (2016)

  * [22] G. Raffelt, D. Dearbom, Bounds on hadronic axions from stellar evolution, Phys. Rev. D36, 2211 (1987)

  * [23] G. Raffelt, Astrophysical methods to constrain axions and other novel particle phenomena, Phys. Rep. 198, 1 (1990)

  * [24] A. Payez, C. Evoli, T. Fischer, M. Giannotti, A. Mirizzi, A. Ringwald, Revisiting the SN1987A gamma-ray limit on ultralight axion-like particles, J. Cosmol. Astropart. Phys. 02, 006 (2015)

  * [25] V. Anastassopoulos et al. (CAST Collaboration), New CAST limit on the axion-photon interaction, Nat. Phys. 13, 584 (2017)