# Faraday effect optical sensing of single-molecules by graphene-based layered
structures

D. Jahani<EMAIL_ADDRESS>O. Akhavan, A. Alidoust Ghatar, W. Fritzsche
and F. Garwe

###### Abstract

_Recently, modulation of the energy bandgap of graphene when gas molecules are
adsorbed to its surface has been proved to be possible. Motivated by this,
based on numerical calculations, we investigate the effect of the associated
bandgap opening in graphene’s spectrum on the sensing properties of the
Faraday rotation (FR) of linearly polarized electromagnetic modes in a 1D
photonic graphene-based sensor with a microcavity defect channel covered by
two graphene layers. Our proposed model introduces an optical contact-free
mechanism for the detection of low-concentration molecules attached to
graphene’s surface. We also show that FR angles could reveal reversal signs
for a different amount of surface transfer doping._

_Department of Physics, Sharif University of Technology, P.O. Box 11155-9161,
Tehran, Iran_

_Leibniz Institute of Photonic Technology (IPHT), Jena, Germany._

_Keywords_ : _Single-molecule detection; Faraday rotation; Graphene; Photonic
Crystal; Fluidic microchannel._

##  1 Introduction

_Graphene, a single sheet of carbon atoms arranged in a honeycomb structure,
was isolated experimentally for the first time from graphite in 2004 [1,2,3].
Thanks to its massless Dirac charge carries, this one-atom-thick material
mimics the properties of a gapless semiconductor and has shown outstanding
electronic and optical properties even at room temperature [4,5,6]. These
interesting properties make graphene serve as an ideal candidate for real-
world technologies. Recently, graphene has also shown huge attention-getting
achievements as an additional key element in many sensing applications
[7,8,9,10]. Surprisingly, because of the low-noise nature of massless charge
carriers of graphene, it an ideal candidate for the detection of individual
gas molecules adsorbed on its surface [11]. Then, the exceptional electronic
property allows the detection of a few molecules on graphene upon surface
transfer doping which correspondingly changes the electrical conductivity of
graphene in order to be measured by proper electrical contacts. In addition to
its unique electrical properties, graphene also reveals peculiar optical
characteristics including, high transparency and broadband ultra-fast response
to the electromagnetic radiation which could interact with its massless Dirac
fermions [12,13]. Consequently, graphene is used in optical contact-free
sensing methods for a broad range of light frequencies [14,15,16,17].
Therefore, optical spectroscopies in the presence of graphene provide a
highly-sensitive nonintrusive method toward the chip-scale sensing technology
in different branches of science. In this regard, in biotechnology, photonic
graphene-based biosensors operating in a relatively low-frequency background
of light prorogation have received considerable attention in recent years
[18,19]. Interestingly, the optical sensing properties of graphene could be
enhanced by applying a constant magnetic field normal to its surface. In fact,
due to induced optical excitations in the surface conductivity of graphene,
applying a magnetic field could rotate the polarization plane of the impinging
light [20]. This finite rotation of the polarization’s plane of the light
propagating through a dielectric medium under Hall situation is called Faraday
rotation (FR) which could be also applied for remote sensing of magnetic field
intensities. Surprisingly, large FR under the Hall regime has been reported
for graphene in spite of being just a one-atom-thick atomic layer [21,22].
This FR in graphene occurs due to intra- and inter-band excitations in its
surface optical conductivity which could be controlled by the magnetic field
intensity [23]. However, there is a trade-off between FR and the
transmission’s height of the transmitted light and, therefore, the possibility
of finding simultaneously both large FR angles and relatively high
transmissions is of high interest in graphene-based integrated optical devices
[24,25,26]._

_Integrated optical technology has provided a broad range of opportunities for
graphene to be used in sensing applications [27,28]. In fact, the presence of
graphene in these devices could improve the sensitivity of optical contact-
free methods which is of prime interest in biotechnology [29,30]. Moreover,
tunability is another important feature for detection purposes since it
suggests more control over the sensing process in optical instruments. In this
regard, 1D photonic crystals (PCs), periodically multilayer structures that
are designed to control the light prorogation could give more tunability
features due to the forbidden frequency zone in the corresponding light
spectrum [31]. Interestingly, creating a defect in these optical structures
could induce a defect mode within the forbidden region which opens a more
controllable route for electromagnetic waves prorogation. Fortunately, growing
successes in photonics applications provide cost-effective components for PCs
in order to host graphene in chemical and biosensing devices. Now, the effect
of the magnetic field on the plane of the polarization of EM waves in a 1D PC
could serve for the detection of gas molecules as dopants on the graphene’s
surface which might show a Dirac gap in its linear dispersion relation. To the
best of our knowledge still little has been reported on FR in defective PCs
for sensing applications in spite of the fact that the Faraday effect itself
is employed for remote sensing of magnetic fields. Moreover, the surface
conductivity of graphene in a magnetic field upon absorption which could also
lead to the creation of Dirac gaps gives different responses to decoupled
polarizations and, therefore, the angle of the rotation could be affected by
the absorption of molecules [32]._

_The first incorporation of graphene into a 1D photonic crystal was introduced
by O. Berman [33,34]. Since then graphene has been offered to be employed in
1D PCs in order to enhance the optical properties of these layer structures.
For instance, the effect of a magnetic field on the polarization plane of
defect states emerged into a 1D graphene-based defective PC proved to enhance
the Faraday effect for which the plane of the polarization of light is rotated
while passing into a magnetic material [35]. Interestingly, the Faraday effect
in graphene-based PCs could be considered in optical sensing devices. In this
work, we address the magneto-sensing properties of a 1D optical structure
composed of a fluidic channel covered by two-graphene layers. In fact, it is
the presence of graphene-coated on internal walls of the fluidic channel which
leads to the rotation of the linearly polarized light considered to be an
equal combination of left- and right-handed circular modes. Note that, to be
applicable in real-world applications, the Faraday effect must be accompanied
by decently high transmission peaks [36]. Also, the realization of a Dirac gap
upon coating graphene and also applying a strong magnetic field normal to its
surface is expected. Therefore, the effect of opening Dirac gaps in graphene’s
spectrum considering the low Fermi energy’s chances which occur due to the
absorbed gas molecules on the graphene’s surface on FR angles is worth huge
attention. To the best of our knowledge, to this stage, no study has been
reported on the effect of opening a Dirac gap on the Faraday effect of
circular modes in a 1D graphene-based defective PC._

_This paper is organized into the following sections: The model and theory is
introduced in section 2. Numerical calculations are addressed in section 3
where we illustrate the effect of the opening an electronic bandgap on the
low-concentration molecules sensing by creation of FR angles in the suggested
optical sensor. Finally, the conclusion is introduced in section 4._

##  2 Model and formalism

_In our model, the Faraday effect of a 1D photonic sensing device which is
schematically shown in figure 1 is investigated based on the transfer matrix
method for a linearly polarized light which is composed of two decoupled
circularly modes impinging normally in $z$-direction. Two graphene sheets in
quantum Hall situation are covering a microfluidic channel, D, in which gas
molecules and other analyzes could be absorbed to graphene’s surface so that
one would expect opening an electronic bandgap in its spectrum which could be
sensed by measuring the corresponding FR angle. The operational principle of
the device relies on the changes in the ultra-sensitive optical conductivity
of graphene due to the absorption of molecules on graphene’s surface.
Interestingly, a short-time ultraviolet illumination could clean the graphene
surface reverting it to its undoped state. Dielectric layers $Si$ and
$SiO_{2}$ with dielectric constants $2.2$ and $10.3$, respectively, form an
optical multilayer system around the graphene-covered microchannel as
$air(AB)^{m}GDG(AB)^{m}air$ for which $m$ indicates the period of the
dielectric layers around the microchannel. The FR angles are then obtained by
calculating the phase difference between transmission coefficients of right-
and left-handed modes. To make things as simple as possible, we consider a
linear polarized light (s or p waves) normally incident on the graphene with
$\varsigma_{s}=1$ and $\varsigma_{p}=-1$:_

$\begin{cases}E_{1s/p}=\varsigma_{s/p}a_{1s/p}e^{ik_{1}z}+b_{1s/p}e^{-ik_{1}z}\\\
H_{1s/p}z_{1}=-a_{1s/p}e^{ik_{1}z}+\varsigma_{s/p}b_{1s/p}e^{-ik_{1}z}\\\
\end{cases}$ (1)

_Now, two circularly polarized modes in the second zone with $E_{2\pm}$ and
$H_{2\pm}$ could be represented as:_

$\begin{cases}E_{2\pm}=a_{2\pm}e^{ik_{2}z}+b_{2\pm}e^{-ik_{2}z}\\\
H_{2\pm}=a_{2\pm}e^{ik_{2}z}+b_{2\pm}e^{-ik_{2}z}\\\ \end{cases}$ (2)

_where $k_{1}=\sqrt{n_{1}}\frac{\omega}{c}$
($k_{2}=\sqrt{n_{2}}\frac{\omega}{c}$) is the associated wave vectors and
$\omega$ and $c$ are the angular frequency and the speed of light in the
vacuum, respectively. Note that $z_{i}$ (i=1,2) is assumed to stand for
$z_{0}/n_{i}$ and $a_{i\pm}=a_{ix}\pm ia_{iy}$ and $b_{i\pm}=b_{ix}\pm
ib_{iy}$. From the boundary condition in the presence of a constant magnetic
field which changes the optical surface conductivity of graphene one might
write:_

$\begin{cases}E_{2\pm}-E_{1s/p}=0\\\ H_{2\pm}-H_{1s/p}=\pm
iJ_{\mp}\end{cases}$ (3)

_for which $J_{\pm}=\sigma_{\mp}E_{2\pm}$. To proceed, we write:_

$\begin{cases}\varsigma_{s/p}a_{1s/p}+b_{1s/p}=a_{2\pm}+b_{2\pm}\\\
-a_{1s/p}+\varsigma_{s/p}b_{1s/p}=(a_{2\pm}+b_{2\pm})-\sigma_{\mp}(a_{2\pm}+b_{2\pm})\\\
\end{cases}$ (4)

_Note that, the longitudinal parts of the optical conductivity tensor $\sigma$
for graphene which proves to be more proper for the numerical calculations are
presented in [37] as:_

$\begin{split}\sigma_{0}(\omega)=&\frac{e^{2}v_{f}^{2}|eB|\left(\hbar\omega+2i\Gamma\right)}{\pi
i}\times\\\
&\sum_{n=0}^{\infty}\left\\{\frac{\left[f_{d}(M_{n})-f_{d}(M_{n+1})\right]+\left[f_{d}(-M_{n+1})-f_{d}(-M_{n})\right]}{\left(M_{n+1}-M_{n}\right)^{3}-\left(\hbar\omega+2i\Gamma\right)^{2}\left(M_{n+1}-M_{n}\right)}\right\\}\\\
&+\left\\{\frac{\left[f_{d}(-M_{n})-f_{d}(M_{n+1})\right]+\left[f_{d}(-M_{n+1})-f_{d}(M_{n})\right]}{\left(M_{n+1}+M_{n}\right)^{3}-\left(\hbar\omega+2i\Gamma\right)^{2}\left(M_{n+1}+M_{n}\right)}\right\\}\end{split}$
(5)

_while the Hall conductivity term is:_

$\begin{split}&\sigma_{H}(\omega)=\frac{-e^{2}v_{f}^{2}eB}{\pi}\\\
&\sum_{n=0}^{\infty}\left\\{\left[f_{d}(M_{n})-f_{d}(M_{n+1})\right]-\left[f_{d}(-M_{n+1})-f_{d}(-M_{n})\right]\right\\}\times\\\
&\left\\{\frac{1}{\left(M_{n+1}-M_{n}\right)^{2}-\left(\hbar\omega+2i\Gamma\right)^{2}}+\frac{1}{\left(M_{n+1}+M_{n}\right)^{2}-\left(\hbar\omega+2i\Gamma\right)^{2}}\right\\}\end{split}$
(6)

_In the above relations $M_{n}=\sqrt{\Delta^{2}+2nv_{F}^{2}|eB|h}$ shows the
quantized energy levels which as it is clear could be tuned by the tunable
Dirac gap. Interestingly, in our case surface transfer doping of graphene by
single-molecules could correspondingly open an electronic bandgap in graphen’s
dispersion. Therefore, these single-molecules could be sensed by the
associated measured FR angles of 1D PC. The distribution are then controlled
by $f_{d}(\varepsilon)=\frac{1}{1+exp(\frac{\varepsilon-\mu_{c}}{K_{B}T)}}$
with $K_{B}$ the Boltzmann constant. The transfer matrix, $D1\rightarrow 2$,
yields the following relation for coupled coefficients of circularly hybrid
waves as:_

$\begin{bmatrix}a_{1\pm}\\\
b_{1\pm}\end{bmatrix}=\frac{1}{2n_{1}}\begin{bmatrix}n_{1}+n_{2}+z_{0}\sigma_{\mp}&&n_{1}-n_{2}+z_{0}\sigma_{\mp}\\\
n_{1}-n_{2}-z_{0}\sigma_{\mp}&&n_{1}+n_{2}-z_{0}\sigma_{\mp}\end{bmatrix}\begin{bmatrix}a_{2\pm}\\\
b_{2\pm}\end{bmatrix}$ (7)

_Therefore, extending this relation to the introduced photonic device, one can
extend it to other layers ( $i-1\rightarrow i$; $i=2,3,4,...$) noticing that
generally complex conductivity $\sigma_{\textpm}$ can be neglected in the
absence of the graphene layer. Hence, for numerical calculations we could use
the following relation:_

$\begin{bmatrix}a_{1\pm}\\\
b_{1\pm}\end{bmatrix}=\frac{1}{2n_{1}}\begin{bmatrix}n_{1}+n_{2}+2\alpha\sigma_{s\mp}&&n_{1}-n_{2}+2\alpha\sigma_{s\mp}\\\
n_{1}-n_{2}-2\alpha\sigma_{s\mp}&&n_{1}+n_{2}-2\alpha\sigma_{s\mp}\end{bmatrix}\begin{bmatrix}a_{2\pm}\\\
b_{2\pm}\end{bmatrix}$ (8)

_Now, in terms of dimensionless terms which are proper for numerical
calculations one can simplify more with the help of the definitions:
$\sigma_{\pm}=\sigma_{0}\sigma_{s\pm}$ and
$\sigma_{0}=\frac{1}{\pi}\frac{e^{2}}{h}$ to show circularly transmission
coefficients for a graphene layer sandwiched between to dielectric media with
refractive indices, $n_{1}$ and $n_{2}$, as:_

$t_{\pm}=\frac{1}{D^{(1,1)}}=\frac{2n_{1}}{\beta+2\alpha\sigma_{s\mp}}$ (9)

_where $\alpha=\frac{e^{2}}{2\pi\epsilon_{0}ch}$ is fine structure constant
and $\beta=n_{1}+n_{2}$. Then, one can show that the transmission coefficients
could be expressed as follows:_

$t_{\pm}=\frac{2n_{1}}{\gamma}[\beta+2\alpha Re(\sigma_{s\mp})-i2\alpha
Im(\sigma_{s\mp})]$ (10)

_with $\gamma=[\beta+2\alpha
Re(\sigma_{s\mp})]^{2}+4\alpha^{2}[Im(\sigma_{s\mp})]^{2}$. Here, one can then
proceed more to obtain FR angle, $\theta_{F}$, in terms of dimentionless
quantities in the following expression:_

$2\theta_{F}=tan^{-1}[-\frac{2\alpha Im(\sigma_{s-})}{\beta+2\alpha
Re(\sigma_{s-})}]-tan^{-1}[-\frac{2\alpha Im(\sigma_{s+})}{\beta+2\alpha
Re(\sigma_{s+})}]$ (11)

_It is clear that for calculating the FR angle of a layered structure we
should also consider a propagation matrix which connects both components of
light (magnetic and electric) to the neighbor dielectric layers:_

$P(z)=\begin{bmatrix}e^{-ik_{z}z}&&0\\\ 0&&e^{ik_{z}z}\end{bmatrix}$ (12)

_At this stage, for the layers structure acting in infrared zone, gathering of
both transmission matrix and propagation one for the first and the last layers
gives :_

$H=M^{(1,2)}P(x^{(1,2)})M^{(2,3)}P(x^{(2,3)})....P(x^{(N-1,N)})M^{(N,N+1)}$
(13)

_The transmission probability for both circularly polarized (left- and right
handed) light are calculated as:_

$T_{\pm}=\frac{1}{H^{(1,1)}}=|t_{\pm}|^{2}$ (14)

_Finally, the transmission probability of the transmitted light is obtained
as:_

$T=\frac{n_{2}}{2n_{1}}(T_{+}+T_{-})$ (15)

_from which FR angle is also calculated based on $\theta_{\pm}$, i.e. the
phase of the left- and right-handed modes as:_

$\theta_{F}=\frac{1}{2}(\theta_{+}-\theta_{-})$ (16)

##  3 Numerical results

_Considering the central wavelength $\lambda_{0}=143\ \mu m$, in this section,
we numerically calculate FR angles and also the transmission probability of
circular microcavity defect modes at room temperature for low values of the
Fermi energy ranging mostly from $\mu=0.1\ meV$ to $\mu=0.9\ meV$ which is
almost $1000$ times smaller than previous works studying the optical
properties of graphene-based sensors. These low Fermi energy changes allow one
to detect low-concentrations of molecules attached to the graphene’s surface.
Moreover, the sign of these changes in the chemical potential could help to
indicate the p-type or n-type doping of attached molecules. We first, as our
main aim in this work, examine the effect of a gap opening on the FR angles of
the linearly polarized defect modes of the suggested optical sensing
structure. As it is revealed in Fig. 2 for $B=10\ T$ and $\mu=0.9\ meV$,
opening a Dirac gap in graphene’s spectrum could significantly affect FR
angles. In fact, increasing the electronic gap at Dirac points in graphene’s
spectrum by attaching low-weight molecules to its surfaces shows to decrease
the rotation of the light’s polarization so that, as it is guessed from the
plots, for $\Delta=200\ meV$ and beyond almost no FR effect is expected to be
observed.
In Fig. 3, the effect of both positive and negative chemical potential on the
FR angles of the defective photonic structure is presented. It is clear that
for a constant magnetic field, $B=10\ T$, and $\Delta=0$ increasing the
chemical potentials gives stronger rotations to the polarization’s plane of
the defect states. Note that, as it is observed, negative values of the
chemical potential cause the linearly polarized light to rotate in an opposite
direction which means that for p-type transfer doping one gets opposite signs
for the rotation angles of the linearly polarized defect modes. Interestingly,
however, we see that the transmission of the defect modes increases as Dirac
gaps open wider and wider. Therefore, as a matter of more illustration, this
behavior is depicted in a 3D plot (see Fig. 4) for opening gaps ranging from
$1$ to $150\ meV$. Interestingly, we see that in the proposed interval for the
electronic bandgap, the transmission probability remains relatively constant.
However, the FR angles exhibit significant changes by increasing the gap
values.
In the next step, in Fig. 5, we simulate the behavior of resonance defect
modes for a constant opening Dirac gap of $50\ meV$ and applied magnetic
fields $B=5$, $B=10$ and $B=20\ T$. Surprisingly, it is seen that stronger
fields could decrease the rotations of the polarization plane of light in some
cases. However, this behavior is not generally true for the transmission of
microcavity resonance modes as it is illustrated. Generally, however, we see
that the effect of increasing the magnetic field on the FR angles is shifting
the angle of rotations toward the central wavelength $\lambda_{0}=143\ \mu m$.
The reason behind shifting the rotation angle toward central frequency modes
is that increasing the intensity of the magnetic field causes separated left-
and right-handed circular modes of the linearly light to move toward merging
in the frequency space.
Furthermore, in Fig. 6, we consider the FR effect for different values of the
chemical potential ranging from $0.1$ to $0.9\ meV$ in the case of gap
openings $\Delta=5$ and $\Delta=50\ meV$. Here, we could see that the
modulation of the chemical potential toward higher values leads to stronger
rotations of the polarization’s plane of the light at a constant magnetic
field $B=10\ T$. It is also obvious that in the case of wider opening for
Dirac gaps in graphene’s spectrum the rotations are observed to be shifted
slightly toward the central wavelength. This situation is better visible fin a
3D view in Fig. 7 for which it is clear that the rotation of the
polarization’s plane of resonance modes shows a linear change relative to the
changes in the value of the chemical potential. It also clear that for higher
values of the chemical potential the rotations of the resonance modes are more
affected by the gap opening.
At this point, we also turn our attention to the case of adding some layers to
our structure. Fig. 8 depicts the situation for the effect of increasing the
gap value on the FR angles in the case of m=7 and a constant magnetic field
$B=10\ T$ at the room temperature. We see that compared to the situations
shown in Fig. 2 we now have stronger rotations. One also could notice the role
of the opening gap in merging the two left- and right-handed modes.
At the end of this section, as a matter of more comparison, we consider the
situation when one of the graphene layers is removed. Interestingly, as is
illustrated in figure 9, both transmission probability and FR angle shown to
be decreased. Particularly, for an opened electronic bandgap in graphene’s
spectrum, $\Delta=60\ meV$ we see that the removal of one of the graphene
layers leads to lower FR angles. On the other hand, as it might have been
expected one gets stronger transmission for the single covering the case of
the microchannel by graphene._

##  4 Conclusion

_To investigate the sensing properties of the Faraday effect in graphene, we
considered the modulations of the energy bandgap of graphene when gas
molecules are adsorbed to its surface. In this work, we numerically
investigated the effect of a tunable electronic gap opening in graphene’s
spectrum on FR rotation for sensing applications of a 1D graphene-based
defective PC. These induced electronic bandgaps could be opened in graphene’s
linear dispersion while single gas molecules are attached to its surface [38].
We demonstrated that the presence of these Dirac gaps affect significantly FR
angles corresponding to the light prorogation into the suggested photonic
sensing device due to the nonlinear effects induced in the energy spectrum of
graphene at Dirac points. On-line monitoring was then proved to be possible by
considering the sign and amount of FR angles in a constant magnetic field
applied on gapped graphene layers for low chemical potentials. To be more
specific, in this work, we first obtained the transmission and the Faraday
effect for the proposed 1D structure with a microfluidic channel covered by
gapped graphene layers in the case of relatively very low changes in the value
of the Fermi energy which could be due to the single molecular absorption on
the graphene’s surface. The existence of an opening gap at Dirac points in
graphene’s spectrum was shown to decrease significantly FR rotations while
keeping the transmission almost unchanged in the proposed photonic sensing
structure. Also, in addition to advances in the electronic bandgap engineering
of graphene, a gap in graphene’s spectrum could be created by applying strong
magnetic fields( excitonic gap) [39]. Therefore, FR angles studied in this
work could provide an optical contact-free method to measure and sense
corresponding opening bandgaps created in the linear graphene’s energy
dispersion relation under different situations. Moreover, it should be pointed
out that all the numerical calculations were carried out at room temperature
and relatively low magnetic fields which could open routes towards the use of
the introduced sensing layered device in real-world applications.
Significantly, we showed that the FR angle for a linearly polarized light
shows narrower bandwidth upon increasing the period of the 1D photonic
structure and also gets relatively higher values by increasing the opening
bandgap. To be more specific, it was observed that inducing an electronic
bandgap $\Delta=50meV$ in graphene’s spectrum for the period $m=7$ gives
stronger FR angles relative to the case for which $m=5$. More importantly,
preservations of Dirac cones upon chemical doping in a very short time
increase the chemical potential of graphene i.e. $\mu=E_{F}-E_{D}$ and thus an
FR occurs for low weight-molecules absorption on graphene in a very short time
could be detected. Therefore, relatively short response times in sensing
applications are expected for the proposed 1D defective photonic structure. In
the end, we showed that the existence of both graphene layers covering the
microchannel cavity is essential since, as we showed numerically, the removal
of one of the graphene layers decreased the FR angles significantly.
In summary, the existence of the bandgap-tunable Faraday effect in a 1D
graphene-based defective optical sensing structure was numerically
demonstrated. Moreover, sensing properties of the introduced model upon
different surface transfer doping - which could lead to opening a Dirac gap in
graphene’s linear spectrum - was studied._

## References

  * [1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Electric field effect in atomically thin carbon films, Science, 306(5696), 666-669.
  * [2] A. C. Neto, F. Guinea, N. M. Peres, K. S. Novoselov and A. K. Geim. The electronic properties of graphene. Reviews of modern physics, 81(1), 109 (2009).
  * [3] A. K. Geim and K. S. Novoselov. The rise of graphene. In Nanoscience and technology: a collection of reviews from nature journals (pp. 11-19) (2010).
  * [4] K. S. Novoselov, Z. Jiang, Y. Zhang, S. V. Morozov, H. L. Stormer, U. Zeitler, J. C. Maan, G. S. Boebinger, P. Kim and A. K. Geim. Room-temperature quantum Hall effect in graphene. Science, 315(5817), 1379-1379 (2007).
  * [5] R. S. Shishir and D. K. Ferry. Intrinsic mobility in graphene. Journal of Physics: Condensed Matter, 21(23), 232204 (2009).
  * [6] L. A Falkovsky. Optical properties of graphene. In Journal of Physics: conference series (Vol. 129, No. 1, p. 012004). IOP Publishing (2008).
  * [7] T. Kuila, S. Bose, P. Khanra, A. K. Mishra, N. H. Kim, J. H. Lee. Recent advances in graphene-based biosensors. Biosensors and bioelectronics, 26(12), 4637-4648 (2011).
  * [8] E. W. Hill, A. Vijayaragahvan and K. Novoselov. Graphene sensors. IEEE Sensors Journal, 11(12), 3161-3170 (2011).
  * [9] K. Z. Milowska and J. A. Majewski. Graphene-based sensors: theoretical study. The Journal of Physical Chemistry C, 118(31), 17395-17401 (2014).
  * [10] E. Vermisoglou, D. Panáček, K. Jayaramulu, M. Pykal, l. Frébort, M. Kolář, M. Hajdúch, R. Zbořila and M. Otyepka. Human virus detection with graphene-based materials. Biosensors and Bioelectronics, 166, 112436 (2020).
  * [11] F. Schedin, A. K. Geim, S. V. Morozov, E. W. Hill, P. Blake, M. I. Katsnelson and K. S. Novoselov. Detection of individual gas molecules adsorbed on graphene. Nature materials, 6(9), 652-655 (2007).
  * [12] D. E. Sheehy and J. Schmalian. Optical transparency of graphene as determined by the fine-structure constant. Physical Review B, 80(19), 193411 (2009).
  * [13] M. Baudisch, A. Marini, J. D. Cox, T. Zhu, F. Silva, S. Teichmann, M. Massicotte, F. Koppens, L. S. Levitov, F. J. García de and Biegert. Ultrafast nonlinear optical response of Dirac fermions in graphene. Nature communications, 9(1), 1-6 (2018).
  * [14] H. Ma, D. Wu, Z. Cui, Y. Li, Y. Zhang, B. Du and Q. Wei. Graphene-based optical and electrochemical biosensors: a review. Analytical Letters, 46(1), 1-17 (2013).
  * [15] M. Hernaez, C. R. Zamarreño, S. Melendi-Espina, L. R. Bird, A. G. Mayes and F. J. Arregui. Optical fibre sensors using graphene-based materials: A review. Sensors, 17(1), 155 (2017).
  * [16] A. A. Rifat, G. A. Mahdiraji, R. Ahmed, D. M. Chow, Y. M. Sua, Y. G. Shee and F. M. Adikan. Copper-graphene-based photonic crystal fiber plasmonic biosensor. IEEE Photonics Journal, 8(1), 1-8 (2015).
  * [17] Z. Li, W. Zhang and F. Xing. Graphene optical biosensors. International journal of molecular sciences, 20(10), 2461 (2019).
  * [18] K. V. Sreekanth, S. Zeng, K. T. Yong and T. Yu. Sensitivity enhanced biosensor using graphene-based one-dimensional photonic crystal. Sensors and Actuators B: Chemical, 182, 424-428 (2013).
  * [19] J. N. Dash and R. Jha. On the performance of graphene-based D-shaped photonic crystal fibre biosensor using surface plasmon resonance. Plasmonics, 10(5), 1123-1131 (2015).
  * [20] I. Fialkovsky and D. V. Vassilevich. Faraday rotation in graphene. The European Physical Journal B, 85(11), 384 (2012).
  * [21] I. Crassee, J. Levallois, A. L. Walter, M. Ostler, A. Bostwick, E. Rotenberg, T. Seyller, D. Van der Marel and A. B. Kuzmenko. Giant Faraday rotation in single-and multilayer graphene. Nature Physics, 7(1), 48-51 (2011).
  * [22] J. C. Martinez, M. B. A. Jalil and S. G. Tan. Optical Faraday rotation with graphene. Journal of Applied Physics, 113(17), 17B529 (2013).
  * [23] V. P. Gusynin, S. G. Sharapov and J. P Carbotte. Magneto-optical conductivity in graphene. Journal of Physics: Condensed Matter, 19(2), 026222 (2006).
  * [24] H. Da, Q. Bao, R. Sanaei, J. Teng, K. P. Loh, F. J. Garcia-Vidal and C. W. Qiu. Monolayer graphene photonic metastructures: Giant Faraday rotation and nearly perfect transmission. Physical Review B, 88(20), 205405 (2013).
  * [25] H. Da and G. Liang. Enhanced Faraday rotation in magnetophotonic crystal infiltrated with graphene. Applied Physics Letters, 98(26), 261915 (2011).
  * [26] H. Da and C. W. Qiu. Graphene-based photonic crystal to steer giant Faraday rotation. Applied Physics Letters, 100(24), 241106 (2012).
  * [27] G., Stijn, G. Navickaite, C. Monasterio, S. Gupta, J. J. Piqueras, R. Pérez, G. Burwell, I. Nikitskiy, T. Lasanta, T. Galán, E. Puma, A. Centeno, A. Pesquera, A. Zurutuza, G. Konstantatos and F. Koppens. Broadband image sensor array based on graphene–CMOS integration. Nature Photonics 11, no. 6 366-371 (2017).
  * [28] X. Gan, C. Zhao, Q. Yuan, L. Fang, Y. Li, J. Yin, X. Ma and J. Zhao. High performance graphene oxide-based humidity sensor integrated on a photonic crystal cavity. Applied Physics Letters 110, no. 15 151107 (2017).
  * [29] Tang, J., Ye, Y., Xu, J., Zheng, Z., Jin, X., Jiang, L. and Xiang, Y. High-Sensitivity Terahertz Refractive Index Sensor in a Multilayered Structure with Graphene. Nanomaterials, 10(3), 500 (2020).
  * [30] Sreekanth, K. V., Zeng, S., Yong, K. T., and Yu, T. Sensitivity enhanced biosensor using graphene-based one-dimensional photonic crystal. Sensors and Actuators B: Chemical, 182, 424-428 (2013).
  * [31] R. D. Meade, J. N. Winn and J. D. Joannopoulos. Photonic crystals: Molding the flow of light. (1995).
  * [32] A. Hosseingholipourasl, S. Hafizah Syed Ariffin, Y. D. Al-Otaibi, E. Akbari, F. Hamid, S. S. R. Koloor and M. Petrů. Analytical Approach to Study Sensing Properties of Graphene Based Gas Sensor. Sensors, 20(5), 1506 (2020).
  * [33] O. L. Berman, V. S. Boyko, R. Y. Kezerashvili, A. A. Kolesnikov and Y. E. Lozovik. Graphene-based photonic crystal. Physics Letters A, 374(47), 4784-4786 (2010).
  * [34] O. L. Berman and R. Y. Kezerashvili. Graphene-based one-dimensional photonic crystal. Journal of Physics: Condensed Matter, 24(1), 015305 (2011).
  * [35] Z. Guo, F. Wu, C. Xue, H. Jiang, Y. Sun, Y. Li, Y and H. Chen. Significant enhancement of magneto-optical effect in one-dimensional photonic crystals with a magnetized epsilon-near-zero defect. Journal of Applied Physics, 124(10), 103104 (2018).
  * [36] H. Da and C. W. Qiu. Graphene-based photonic crystal to steer giant Faraday rotation. Applied Physics Letters, 100(24), 241106 (2012).
  * [37] V. P. Gusynin, S. G. Sharapov and J. P. Carbotte. Magneto-optical conductivity in graphene. Journal of Physics: Condensed Matter, 19(2), 026222 (2006).
  * [38] A. Hosseingholipourasl, S. Hafizah Syed Ariffin, Y. D. Al-Otaibi, E. Akbari, F. Hamid, S. S. R. Koloor and M. Petrů. Analytical Approach to Study Sensing Properties of Graphene Based Gas Sensor. Sensors, 20(5), 1506 (2020).
  * [39] V. P. Gusynin, V. A. Miransky, S. G. Sharapov and L. A. Shovkovy . Excitonic gap, phase transition, and quantum Hall effect in graphene. Physical Review B, 74(19), 195429 (2006).

Figure 1:  The proposed defective optical sensing device with a fluidic
microchannel covered by graphene layers. The changes in FR angles could be
observed by attaching low weight molecules to gapped-graphene’s surface.
Figure 2: Effect of gap opening at Dirac points in graphene on Faraday modes
and the corresponding transmission of EM waves. It is obvious that there is a
trade-off between the values of electronic gaps and FR effect in the
introduced optical sensing structure with the period number $m=5$.  Figure 3:
Faraday modes for gapless graphene are shown to posses both positive and
negative signs which corresponds to the two direction of FR angles for the
different doping transfer on the graphene’s surface due to changes in the sign
of the chemical potential. It also observed that increasing the chemical
potential leads to the observation of stronger rotations of the polarization’s
plane of EM modes. Figure 4: Effect of gap opening $\Delta$ at Driac points in
graphene’s energy dispersion on FR angles and also the transmission modes for
a constant chemical potential $\mu=0.9\ meV$. Interestingly, as it is seen, in
the interval $\Delta=0$ to $\Delta=150\ meV$ for the electronic bandgap, FR
angles show to be more sensitive to the changes of $\Delta$ than the
transmission probability which is relatively constant. Figure 5: The frequency
shift of FR angles for two different Dirac gaps is illustrated. It is also
shown that the effect of increasing the chemical potential for the both
situation leads to no frequency shift for FR angles. Figure 6: The effect of
increasing chemical potential at $B=10\ T$ for gaps $\Delta=5$ and $\Delta=50\
meV$. The blue shift of the modes are observed for the higher Dirac opening.
Figure 7: Faraday effect for different values of Dirac gap opening at room
temperature. Linear relation between FR angles and the value of the chemical
potential is observed in this 3D plot for two different $\Delta=5\ meV$ and
$\Delta=50\ meV$. Significantly, it is seen that the separation of FR modes
for the two values occurs at higher chemical potentials.  Figure 8: Effect of
gap opening for Fr angles and their corresponding transmissions at the room
temperature for an applied magnetic field $B=10\ T$. As it seen Faraday effect
is more diminished upon higher bandgaps in graphene’s spectrum. Figure 9: For
an opened electronic bandgap, $\Delta=60\ meV$ we see that the removal of one
graphene layer leads to lower FR angles. On the other hand, as it might have
been expected one gets stronger transmission for the single covering case of
the microchannel by graphene.


# Thermodynamic Stable Site for Interstitial alloy (N or O) in bcc-Refractory
Metals using Density Functional Theory

Henry Martin  Peter Amoako-Yirenkyi
hmartin<EMAIL_ADDRESS>Eric K. K. Abavare  Center
for Scientific and Technical Computing, National Institute for Mathematical
Sciences, Kumasi, Ghana  Department of Physics, Kwame Nkrumah University of
Science and Technology, Kumasi, Ghana  Department of Mathematics, Kwame
Nkrumah University of Science and Technology, Kumasi, Ghana

###### Abstract

Plasticity in body centered cubic (bcc) refractory metals are largely due to
the stress tensor induced either by solute or thermal activation. The
mechanism of the solute atom(s) residence causes instability in such metals.
Earlier research have considered the mechanism of oxygen (O) or carbon (C) in
tungsten (W), even though the major component of the environment is nitrogen
(N). In this article, the density functional theory (DFT) was employed to
investigate the thermodynamic stable site for an interstitial solute (N or O)
in the bcc refractory metals (Mo and Nb) by calculating the equilibrium and
structural parameters, dissolution energetics and volumetric strain. The
dissolution mechanism of all the relaxed solid solution structures were
predicted to be an exothermic reaction from the supersaturated cell to the low
concentration (1.82 at.%) except for Mo-N solid solution. Convergence of
volumetric strain was observed at the low concentration of the solute. At this
point, the solid solution of Mo-N and Mo-O had a less measure of global stress
(less distortion) at the octahedral (o) site while that of Nb-N and Nb-O were
at the tetrahedral (t) site. This certainly shows why this two bcc-refractory
metals in groups VB (Nb) and VIB (Mo) of the periodic table exhibit different
deformation behaviours giving their difference in site preference stability.

###### keywords:

Stable site , Thermodynamics , Solid solution , bcc-refractory-interstitial ,
DFT

††journal: Physica B: Condensed Matter

##  1 Introduction

Body centered cubic (bcc) refractory metals have extraordinary resistance to
wear and heat giving their unusual high melting point [1]. They also have high
hardness value at room temperature and are chemically inert with relatively
high density [2]. These metals are relatively abundant naturally in the
earth’s crust as mineral compounds of oxides and sulfides, Figure 1 depicts
few of such elements [3]. The Niobium mineral: _Pyrochlore_
($(Ca,Na)_{2-m}Nb_{2}O_{6}(O,OH,F)_{1-n}xH_{2}O$) and _Columbite_
($(Fe,Mn)(Nb,Ta)_{2}0_{6}$) are forms for obtaining its ore deposits. The only
sulfide compound mineral is _Molybdenite_ ($MoS_{2}$), which is found as veins
in quartz rock [3].

Figure 1: Natural forms of some refractory metals with relative percentage
abundance aligned in groups VB and VIB of the periodic table [3]

According to the Kepler conjecture, bcc crystal structures are less dense in
comparison with fcc (ccp) and hcp [4, 5]. This, at low temperature warrants
different behaviours in these structures during plastic deformation [6, 7, 8].
Particularly at decreasing temperature, there exist rapid increase in the flow
stress and strain rate of single crystals (thus, bcc metals) [9, 10, 11, 12]
which is mainly caused by breakdown in the standard geometric projection (non-
glide components added to the stress tensor) [13, 14] during plastic
deformation by solute [15, 16] or thermal activation [17, 18, 19].
Many researchers have established that the bcc refractory metals of groups VB
(V, Nb, and Ta) and VIB (Cr, Mo, and W) exhibit different deformation
behaviours [20, 21, 22]. This difference is ascribed to the different
behaviours of domination of screw dislocation motion at the cores during
plastic deformation at low temperatures [20, 23, 24]. Woodward and Rao
demonstrated that the periodic table elements of group VB (Nb) has screw
displacement spreading evenly about a central point while group VIB (Mo) are
off centred [25, 26, 27, 28].
At lower temperatures, these bcc refractory metals display poor fracture
toughness. Hence, the need to consider alloying with other elements to improve
their ductility. A known technique called solid solution for producing novel
materials has played significant role from the Bronze Age through the Iron Age
to these advanced technological times [29, 30]. This technique has so far
caused a significant change in the application of bcc metals strengthening
(softening or hardening). [31, 32, 33]. From the onset of this solid solution
technique, several researches (experimentally and computationally) have sought
to investigate and understand substitutional solid solution [16, 34, 17, 35,
18, 19, 36, 37, 38, 39, 40]. Unfortunately, little attention has been given to
the interstitial solid solution [41, 42, 43].
However, some research has demonstrated that interstitial solutes such as B,
C, N, O in bcc transition metal (Fe) stabilizes the hard core configuration of
screw dislocation [38, 28, 44, 45]. Recently, similar research has been done
in bcc-refractory metal for C in (V, Nb, Ta, Mo and W) and O in W [41, 43].
These interstitial solutes have also been proven to produce a marked increase
in the ductile-to-brittle transition temperature (DBTT) of bcc metals. This
relation to DBTT effect is linked to the dislocation motion at a specific
temperature level. A reduction of the kinetic energy (cause of low
temperature) of the atoms slows the movement of dislocation in order for
plasticity to take place. Typically, the addition of O and N to high pure
crystals of Nb and Ta causes hardening [46, 47]. However, research show that
the hardness of Nb decrease with increasing Nb concentration [48, 49, 50, 51].
This shows the presence of the interstitial solutes retards the dislocation
motion for the occurrence of plasticity.
Given the importance of these bcc refractory interstitial solid solution in
various applications where their plastic response ranges across length, time,
and temperature scales. It is critical that we understand their performance
and stability in any applicable regime on any scale [52, 53].
This paper addresses the thermodynamic mechanism occupancy preference of the
interstitial solutes in bcc refractory metals. This was accomplished by using
DFT because of the experimental difficulty of studying and measuring the
effects of atomic distributions and external stresses on thermodynamic
properties [54]. The use of DFT offers an opportunity such as accessing the
compositional temperature range of interest. This grants the ability to
describe the deformation as well as to provide the essential chemical details
of the interaction. Characteristics such as the equilibrium geometries,
dissolution energetics and volumetric strain of the systems (MoN, NbN, NbO and
MoO) were the center of focus for our investigation and analysis.

##  2 Computational Method

The Density Functional Theory (DFT) formulation [55, 56] as implemented in the
Vienna ab initio Simulation Package (VASP) [57, 58, 59] was employed. In
employing DFT in the VASP package a construction of different sizes of
supercell using a range of $1\times 1\times 1$ to $3\times 3\times 3$ was
first done to depict different concentration of solute in the bcc refractory
metal. This was done such that one interstitial solute (N or O) is inserted
within an octahedral or tetrahedral sites of the bulk system (supercell of Mo
or Nb) as shown in Figure 2.

Figure 2: Octahedral and tetrahedral sites of bcc unit cell [43]

The Projector Augmented Wave (PAW) pseudopotential scheme [60, 61] with semi-
core p states was considered as part of the valence electrons for the bcc
refractory metals [62, 63, 64]. The Perdew-Burke-Ernzerhof (PBE) functional
[65, 66] within the Generalized Gradient Approximation (GGA) [67], together
with the Hermite-Gaussian within the Methfessel and Paxton method [68, 69, 70]
were employed with an appropriate smearing sigma to obtain a negligible energy
difference value less than 1 meV$/$atom.
The convergence of kinetic energy cut-off was found with appropriate k-points
leading to a stable thermodynamic system. For an accurate and validated
comparison, the kinetic energy cut-off and the k-point chosen for a particular
supercell structure (example Mo) is used for all systems except that of the
molecules ($\ch{A2}=\ch{N2}$ or $\ch{O2}$). The systems used k-point sampling
of $7\times 7\times 7$, $3\times 3\times 3$ and $3\times 3\times 3$ for
$X_{2}$, $X_{16}$ and $X_{54}$ respectively. Where, $X$ represent all systems.
The molecules were computed by placing the dimer in a cubic box after which a
$\Gamma$ point calculation was carried out. The modified Broyden’s method was
applied for the final step of charge mixing scheme [71]. The relaxation of
both cell shape and atomic positions structures were performed until forces on
each atom was less than $10^{-5}$ eV$/$atom using the Conjugate Gradient (CG)
technique.
The total energies of the interstitial solutes and the bulk refractory metals
were obtained to calculate the dissolution energy or the heat of solution at a
specific x-site ($x=\\{o$ or $t\\}$), $E_{d_{x}}$ as

$E_{d_{x}}=E_{\text{bcc cell}+\text{IS}}-E_{\text{bcc
cell}}-\dfrac{1}{2}E[\ch{A2}]$ (1)

Where $E_{\text{bcc cell}+\text{IS}}$ and $E_{\text{bcc cell}}$ are the
energies for bulk bcc refractory metal with and without an inserted
interstitial solute respectively. The last term defines the energy of the
molecular atom being considered. The dissolution energy difference, $\Delta
E_{d_{o-t}}$, with respect to the difference in the occupied site ($x=\\{o$ or
$t\\}$) was also obtained. Finally, the volume of the supercells were obtained
after the relaxation of the whole system with and without the interstitial
solute. This was done to calculate the volumetric strain which is equivalent
to hydrostatic strain as follows:

$\epsilon_{v}=\dfrac{\Delta V}{V}\approx\operatorname{Tr}(\epsilon)=3\epsilon$
(2)

##  3 Results and Discussion

In determining the thermodynamic stable site of bcc refractory interstitial
solid solution using the DFT technique, the equilibrium geometries were first
obtained. The lattice constant of the bulk bcc refractory metals (Mo and Nb)
and the vibration frequency of the dimers ($N_{2}$ and $O_{2}$) are shown in
Table 1. These as shown in the same Tables are within the range of
experimental data as well as previously calculated DFT results.

Table 1: PAW-DFT-GGA predictions of lattice parameters (a) of bcc refractory
metals (Mo and Nb) and the vibration frequency ($cm^{-1}$) of the dimers (N2
and O2 gas) together with some other computational and experimental results
are also shown for comparison

Specie | a (Å) | Specie | a (Å) | Dimers |  | Freq ($cm^{-1}$) |
---|---|---|---|---|---|---|---
Mo | $3.0996^{a}$ | Nb | $3.294^{g}$ | $N_{2\ {\tiny(gas)}}$ |  | $2359^{m}$ |
$3.142\pm 0.03^{b}$ | $3.296^{h}$ |  | $2566.28^{n}$ |
$3.15^{c,d}$ | $3.300^{i,j}$ |  | $\mathbf{2579.48}^{q}$ |
$3.16^{e}$ | $3.312^{k}$ | $O_{2\ {\tiny(gas)}}$ |  | $1556^{o}$ |
$3.17^{f}$ | $3.326^{l}$ |  | $\mathbf{1645.1}^{q}$ |
$\mathbf{3.17}^{q}$ | $\mathbf{3.327}^{q}$ |  | $2359^{p}$ |

  * $a$

Ref.[16]

  * $b$

Ref.[72]

  * $c$

Ref.[73]

  * $d$

Ref.[74]

  * $e$

Ref.[75]

  * $f$

Ref.[63]

  * $g$

Ref.[73]

  * $h$

Ref.[76]

  * $i$

Ref.[77]

  * $j$

Ref.[78]

  * $k$

Ref.[79]

  * $l$

Ref.[80]

  * $m$

Ref.[81]

  * $n$

Ref.[82]

  * $o$

Ref.[83]

  * $p$

Ref.[84]

  * $q$

This work

As displayed in Figure 2, the experimental lattice parameter ($a$) and atomic
radius ($r$) of the bcc refractory metals (Mo is $3.15$ Å and $1.363$ Å and Nb
is $3.30$ Å and $1.430$ Å respectively) were replaced in a theoretical
formulae given as:

$\frac{r_{o}}{r}=\dfrac{a(\frac{1}{2}-\frac{\sqrt{3}}{4})}{a\frac{\sqrt{3}}{4}}$
(o-site)
$\dfrac{r_{t}}{r}=\dfrac{a(\frac{\sqrt{5}}{4}-\frac{\sqrt{3}}{4})}{a\frac{\sqrt{3}}{4}}$
(t-site)

to obtain the radii of the o-site ($r_{o}$) and t-site ($r_{t}$) for Mo as
$0.21086$ Å and $0.39663$ Å, respectively and for Nb as $0.22122$ Å and
$0.41612$ Å, respectively. N and O has a covalent radius of $0.74$ Å and
$0.66$ Å respectively. Therefore, the expectation would be that N and O will
fit better in the t-site with less distortion, given that the radius at the
t-site is wider than that of the o-site for fixed atoms position or Unrelaxed
structure. We note that this theoretical expectation do not account for any
form of relaxation within the geometry of the bcc lattice.
In confirming this assertion from the theoretical formulae, we first assess
Table 2 \- 5 showing the dissolution energy difference of N and O occupancy in
the o-site and t-site of bcc refractory metals supercell sizes using different
concentration of N and O. Since the difference in dissolution energy
determines the minimum energy at a particular site, the stability is gained
for the alloying element located in the pure metal (bcc-refractory
interstitial solid solution). It can be expressed as:

$\Delta E_{d_{o-t}}=+\ /\ -$ (3)

The ($+$) energy value $\Rightarrow$ $E_{d_{o}}>E_{d_{t}}$; $\therefore$
t-site has minimum energy while the ($-$) energy value $\Rightarrow$
$E_{d_{o}}<E_{d_{t}}$; $\therefore$ o-site has minimum energy. This confirms
the results of the theoretical formulae by considering the unrelaxed
dissolution energy difference of all the bcc-refractory interstitial solid
solution from Table 2 \- 5 as well as the dissolution energies at the Appdx.
A1 \- A4. Unfortunately, the relaxed structure was not permitted to conform to
the same circumstance under which the theoretical formulae was set up. This is
clearly shown for all concentration of the supercells or atleast either the
1st and 2nd supercells for the solid solutions.

Table 2: Dissolution energy difference ($\Delta E_{d_{o-t}}$) between Nitrogen (N) in the o-site and t-site of bcc Molybdenum (Mo) for both unrelaxed (Ur) and relaxed (R) structures Supercell | $C_{N}(at.\%)$ |  $\Delta E_{{Ur}_{o-t}}$ (eV) |  $\Delta E_{R_{o-t}}$ (eV)
---|---|---|---
Mo2N1 | 33.3 | 2.0274 | -1.4906
Mo16N1 | 5.88 | 1.8335 | -0.5429
Mo54N1 | 1.82 | 2.0287 | -0.6960

Starting with the concentration of $X_{2}A$ cells which may be considered as a
supersaturated solution of the interstitial in the bcc refractory metals. Its
noticed that there was very high change in the dissolution energy towards the
negative as compared to the other concentrations. This depicts stability at
the octahedral sites which is in agreement with Table 6 and Figure 3 showing
less volumetric strain at the o-site than that of the t-sites.

Table 3: Dissolution energy difference ($\Delta E_{d_{o-t}}$) between Nitrogen (N) in the o-site and t-site of bcc Niobium (Nb) for both unrelaxed (Ur) and relaxed (R) structures Supercell | $C_{N}(at.\%)$ |  $\Delta E_{{Ur}_{o-t}}$ (eV) |  $\Delta E_{R_{o-t}}$ (eV)
---|---|---|---
Nb2N1 | 33.3 | 1.3852 | -1.6695
Nb16N1 | 5.88 | 1.2986 | -0.9776
Nb54N1 | 1.82 | 1.4313 | 0.0041
Table 4: Dissolution energy difference ($\Delta E_{d_{o-t}}$) between Oxygen (O) in the o-site and t-site of bcc Molybdenum (Mo) for both unrelaxed (Ur) and relaxed (R) structures Supercell | $C_{O}(at.\%)$ |  $\Delta E_{{Ur}_{o-t}}$ (eV) |  $\Delta E_{R_{o-t}}$ (eV)
---|---|---|---
Mo2O1 | 33.3 | 1.7413 | -0.5892
Mo16O1 | 5.88 | 2.7161 | 0.3629
Mo54O1 | 1.82 | 2.8144 | 0.0893
Table 5: Dissolution energy difference ($\Delta E_{d_{o-t}}$) between Oxygen (O) in the o-site and t-site of bcc Niobium (Nb) for both unrelaxed (Ur) and relaxed (R) structures Supercell | $C_{O}(at.\%)$ |  $\Delta E_{{Ur}_{o-t}}$ (eV) |  $\Delta E_{R_{o-t}}$ (eV)
---|---|---|---
Nb2O1 | 33.3 | 1.4803 | -0.1599
Nb16O1 | 5.88 | 1.4312 | -0.8473
Nb54O1 | 1.82 | 1.6656 | 0.0336

Unlike the unrelaxed structures, there is a change with the preference of the
site occupancy when the concentration diminishes. This change of site
preference from the o-site to the t-site in the concentration of $X_{16}A$
supercell is ONLY recognised in Mo16O1 (Table 4) which is clearly seen among
the other bcc refractory interstitial solid solution. It also shows a clear
agreement with the volumetric strain shown in Table 6 even though the Mo16O1
depicts otherwise. This dispute with the dissolution energy difference and the
volumetric strain, may be due to the slight change in shift or distortion as
seen in Table 6, which show that the volumetric strain at the o-site and
t-site are quite closer as compared to other solid solution.

Table 6: Volumetric strain, $\epsilon_{v}=(\frac{\Delta V}{V})^{*}$ for the
relaxed (R) cells structure of the interstitial solute, ($A=\\{N$ or $O\\}$)
in the Octahedral (o) or Tetrahedral (t) site of the bcc refractory metal
($X=\\{Mo$ or $Nb\\}$)

$C_{A}(at.\%)$ | ${Mo}_{x}N1$ | ${Nb}_{x}N1$ | ${Mo}_{x}O1$ | ${Nb}_{x}O1$
---|---|---|---|---
$\epsilon_{v_{o}}$ | $\epsilon_{v_{t}}$ | $\epsilon_{v_{o}}$ | $\epsilon_{v_{t}}$ | $\epsilon_{v_{o}}$ | $\epsilon_{v_{t}}$ | $\epsilon_{v_{o}}$ | $\epsilon_{v_{t}}$
33.3 | 17.0426 | 28.0075 | 13.1502 | 25.7466 | 18.8283 | 25.3133 | 13.5845 | 30.0358
5.88 | 3.1029 | 3.9317 | 2.7804 | 4.2476 | 3.6135 | 3.9709 | 2.9754 | 3.7825
1.82 | 0.9074 | 1.0610 | 1.0028 | 0.8207 | 1.0063 | 1.1168 | 1.0160 | 0.6871

  * $*$

Volume increase of the relaxed $X_{x}A$ cell with respect to the relaxed
$X_{x}$ cell.

The last supercell shows a complete change to the original site preference
proposed by the theoretical formulae (thus, Eqs. o-site and t-site) as
discussed earlier for all the solid solution except for Mo-N. This complete
change again can be linked to the volumetric strain which is in total
agreement with the assertion by taking a look at Figure 3 but has a slight
contrast as shown in Table 6 with respect to the solid solution Mo-O. This
slight contrast between the dissolution energy difference and the volumetric
strain, is due to the closeness of the distortion that occur at both sites
within the Mo54O1. Figure 3 shows the form of convergence with respect to the
distortion and stability when the site preferences are of a concern to all
solid solution.

Figure 3: A plot of volumetric strains when an interstitial solute (N or O) is
inserted in the Octahedral (o) or Tetrahedral (t) site of a bcc refractory
metal (Mo or Nb) at 0 K and relaxed

##  4 Conclusion

Density Functional Theory calculation was performed to examine the
thermodynamic stable site preference on some selected bcc refractory
interstitial solid solution. The dissolution mechanism of all the relaxed
solid solution structures were predicted to be an exothermic reaction from the
supersaturated cell to the low concentration (1.82 at.%) except for Mo-N solid
solution. Convergence of volumetric strain was observed moving towards low
concentration (1.82 at.%) of the solute with a less measure of global stress
(less distortion) at the octahedral (o) site for the solid solutions of Mo-N
and Mo-O and at the tetrahedral site for the Nb-N and Nb-O. This certainly
shows why the two (2) bcc-refractory metals in groups VB (Nb) and VIB (Mo) of
the periodic table exhibit different deformation behaviours giving their
difference in site preference stability. We also noted that the demonstration
of stability for all the solid solution at less concentration requires further
study on the movement of the solute (N or O) across the sites (thus through
o-sites, t-sites and o-t-sites) with a very low concentration of the solute or
bigger supercell than what was under study in this work.

##  5 Acknowledgements

Centre for Scientific and Technical Computing, National Institute for
Mathematical Sciences hosted by Kwame Nkrumah University of Science and
Technology Kumasi, Ghana. The authors would like to thank Prof. Jaime Marian
and his group at Department of Materials Science and Engineering and
Mechanical and Aerospace Engineering, University of California Los Angeles
(UCLA), Los Angeles, USA, for his helpful discussion, insightful advice,
access and usage of the computational and storage services associated with the
Hoffman2 Shared Cluster provided by UCLA Institute for Digital Research and
Education’s Research Technology Group.

##  6 Appendix

Table A1: Dissolution energy ($E_{d}$) of Nitrogen (N) in the Octahedral (o) and Tetrahedral (t) sites of bcc Molybdenum (Mo) for both unrelaxed (Ur) and relaxed (R) structures Supercell | $C_{N}(at.\%)$ | o site | t site
---|---|---|---
$E_{d_{Ur}}$ (eV) |  $E_{d_{R}}$ (eV) |  $E_{d_{Ur}}$ (eV) |  $E_{d_{R}}$ (eV)
Mo2N1 | 33.3 | 2.9713 | -0.6064 | 0.9439 | 0.8842
Mo16N1 | 5.88 | 5.4796 | 0.9791 | 3.6461 | 1.5220
Mo54N1 | 1.82 | 6.0527 | 0.7795 | 4.0240 | 1.4755
Table A2: Dissolution energy ($E_{d}$) of Nitrogen (N) in the Octahedral (o) and Tetrahedral (t) sites of bcc Niobium (Nb) for both unrelaxed (Ur) and relaxed (R) structures Supercell | $C_{N}(at.\%)$ | o site | t site
---|---|---|---
$E_{d_{Ur}}$ (eV) |  $E_{d_{R}}$ (eV) |  $E_{d_{Ur}}$ (eV) |  $E_{d_{R}}$ (eV)
Nb2N1 | 33.3 | 1.0690 | -2.1949 | -0.3162 | -0.5253
Nb16N1 | 5.88 | 2.0490 | -1.7504 | 0.7504 | -0.7728
Nb54N1 | 1.82 | 2.2927 | -2.2332 | 0.8614 | -2.2373
Table A3: Dissolution energy ($E_{d}$) of Oxygen (O) in the Octahedral (o) and Tetrahedral (t) sites of bcc Molybdenum (Mo) for both unrelaxed (Ur) and relaxed (R) structures Supercell | $C_{N}(at.\%)$ | o site | t site
---|---|---|---
$E_{d_{Ur}}$ (eV) |  $E_{d_{R}}$ (eV) |  $E_{d_{Ur}}$ (eV) |  $E_{d_{R}}$ (eV)
Mo2O1 | 33.3 | 1.0479 | -1.4403 | -0.6934 | -0.8512
Mo16O1 | 5.88 | 4.7204 | -0.0226 | 2.0043 | -0.3855
Mo54O1 | 1.82 | 5.3028 | -0.3552 | 2.4885 | -0.4444
Table A4: Dissolution energy ($E_{d}$) of Oxygen (O) in the Octahedral (o) and
Tetrahedral (t) sites of bcc Niobium (Nb) for both unrelaxed (Ur) and relaxed
(R) structures

Supercell | $C_{N}(at.\%)$ | o site | t site
---|---|---|---
$E_{d_{Ur}}$ (eV) |  $E_{d_{R}}$ (eV) |  $E_{d_{Ur}}$ (eV) |  $E_{d_{R}}$ (eV)
Nb2O1 | 33.3 | -1.1301 | -3.5987 | -2.6103 | -3.4388
Nb16O1 | 5.88 | 0.3354 | -3.6867 | -1.0958 | -2.8394
Nb54O1 | 1.82 | 0.6030 | -4.0744 | -1.0627 | -4.1080

  * 1.

($+$) pos energy value - endothermic process of dissolution.

  * 2.

($-$) neg energy value - exothermic process of dissolution

## References

## References

  * [1] J. R. Davis, Alloying: understanding the basics, ASM international, 2001.
  * [2] M. Bauccio, ASM metals reference book ”Refractory metals”, ASM international, 1993\.
  * [3] F. Habashi, Historical introduction to refractory metals, Mineral Procesing and Extractive Metallurgy Review 22 (1) (2001) 25–53.
  * [4] J. Kepler, C. G. Hardie, B. J. Mason, L. L. Whyte, The Six-cornered Snowflake Lat. & Eng, Clarendon Press, 1966.
  * [5] T. Hales, M. Adams, G. Bauer, T. D. Dang, J. Harrison, H. Le Truong, C. Kaliszyk, V. Magron, S. McLaughlin, T. T. Nguyen, et al., A formal proof of the kepler conjecture, in: Forum of Mathematics, Pi, Vol. 5, Cambridge University Press, 2017.  doi:10.1017/fmp.2017.1.
  * [6] J. Christian, Some surprising features of the plastic deformation of body-centered cubic metals and alloys, Metallurgical Transactions A 14 (7) (1983) 1237–1256.  doi:10.1007/BF02664806.
  * [7] D. Caillard, J.-L. Martin, Thermally activated mechanisms in crystal plasticity, Vol. 8, Elsevier, 2003.
  * [8] A. S. Argon, Strengthening mechanisms in crystal plasticity, Vol. 4, Oxford University Press, 2008.
  * [9] A. Seeger, The flow stress of high-purity refractory body-centred cubic metals and its modification by atomic defects, Le Journal de Physique IV 5 (C7) (1995) C745–C765.  doi:10.1051/jp4:1995704.
  * [10] W. Pichl, Slip geometry and plastic anisotropy of body-centered cubic metals, physica status solidi (a) 189 (1) (2002) 5–25.  doi:10.1002/1521-396X(200201)189:1<5::AID-PSSA5>3.0.CO;2-D.
  * [11] A. Stukowski, D. Cereceda, T. D. Swinburne, J. Marian, Thermally-activated non-schmid glide of screw dislocations in w using atomistically-informed kinetic monte carlo simulations, International Journal of Plasticity 65 (2015) 108–130.  doi:10.1016/j.ijplas.2014.08.015.
  * [12] D. Cereceda, M. Diehl, F. Roters, D. Raabe, J. M. Perlado, J. Marian, Unraveling the temperature dependence of the yield strength in single-crystal tungsten using atomistically-informed crystal plasticity calculations, International Journal of Plasticity 78 (2016) 242–265.  doi:10.1016/j.ijplas.2015.09.002.
  * [13] J. Chaussidon, M. Fivel, D. Rodney, The glide of screw dislocations in bcc fe: Atomistic static and dynamic simulations, Acta materialia 54 (13) (2006) 3407–3416.  doi:10.1016/j.actamat.2006.03.044.
  * [14] S. Brinckmann, J.-Y. Kim, J. R. Greer, Fundamental differences in mechanical behavior between two types of crystals at the nanoscale, Physical review letters 100 (15) (2008) 155502.  doi:10.1103/PhysRevLett.100.155502.
  * [15] P. Penning, Mathematics of the portevin-le chatelier effect, Acta Metallurgica 20 (10) (1972) 1169–1175.  doi:10.1016/0001-6160(72)90165-4.
  * [16] D. R. Trinkle, C. Woodward, The chemistry of deformation: How solutes soften pure metals, Science 310 (5754) (2005) 1665–1667.  doi:10.1126/science.1118616.
  * [17] L. Romaner, C. Ambrosch-Draxl, R. Pippan, Effect of rhenium on the dislocation core structure in tungsten, Physical review letters 104 (19) (2010) 195503.  doi:10.1103/PhysRevLett.104.195503.
  * [18] H. Li, S. Wurster, C. Motz, L. Romaner, C. Ambrosch-Draxl, R. Pippan, Dislocation-core symmetry and slip planes in tungsten alloys: Ab initio calculations and microcantilever bending experiments, Acta Materialia 60 (2) (2012) 748–758.  doi:10.1016/j.actamat.2011.10.031.
  * [19] G. D. Samolyuk, Y. Osetsky, R. Stoller, The influence of transition metal solutes on the dislocation core structure and values of the peierls stress and barrier in tungsten, Journal of Physics: Condensed Matter 25 (2) (2012) 025403\.  doi:10.1088/0953-8984/25/2/025403.
  * [20] M. S. Duesbery, V. Vitek, Plastic anisotropy in bcc transition metals, Acta Materialia 46 (5) (1998) 1481–1492.  doi:10.1016/S1359-6454(97)00367-4.
  * [21] C. R. Weinberger, B. L. Boyce, C. C. Battaile, Slip planes in bcc transition metals, International Materials Reviews 58 (5) (2013) 296–314.  doi:10.1179/1743280412Y.0000000015.
  * [22] C. R. Weinberger, G. J. Tucker, S. M. Foiles, Peierls potential of screw dislocations in bcc transition metals: Predictions from density functional theory, Physical Review B 87 (5) (2013) 054114.  doi:10.1103/PhysRevB.87.054114.
  * [23] H. Lim, C. R. Weinberger, C. C. Battaile, T. E. Buchheit, Application of generalized non-schmid yield law to low-temperature plasticity in bcc transition metals, Modelling and Simulation in Materials Science and Engineering 21 (4) (2013) 045015.  doi:10.1088/0965-0393/21/4/045015.
  * [24] H. Lim, C. C. Battaile, C. R. Weinberger, Simulating dislocation plasticity in bcc metals by integrating fundamental concepts with macroscale models, Integrated Computational Materials Engineering (ICME) for Metals: Concepts and Case Studies (2018) 71–105doi:10.1002/9781119018377.ch4.
  * [25] C. Woodward, S. I. Rao, Ab-initio simulation of isolated screw dislocations in bcc mo and ta, Philosophical Magazine A 81 (5) (2001) 1305–1316.  doi:10.1080/01418610108214442.
  * [26] L. Dezerald, L. Proville, L. Ventelon, F. Willaime, D. Rodney, First-principles prediction of kink-pair activation enthalpy on screw dislocations in bcc transition metals: V, nb, ta, mo, w, and fe, Physical Review B 91 (9) (2015) 094105\.  doi:10.1103/PhysRevB.91.094105.
  * [27] L. Dezerald, D. Rodney, E. Clouet, L. Ventelon, F. Willaime, Plastic anisotropy and dislocation trajectory in bcc metals, Nature communications 7 (2016) 11695\.  doi:10.1038/ncomms11695.
  * [28] D. Rodney, L. Ventelon, E. Clouet, L. Pizzagalli, F. Willaime, Ab initio modeling of dislocation core properties in metals and semiconductors, Acta Materialia 124 (2017) 633–659.  doi:10.1016/j.actamat.2016.09.049.
  * [29] M. Hansen, K. Anderko, H. Salzberg, Constitution of binary alloys, Journal of the Electrochemical Society 105 (12) (1958) 260C–261C.  doi:10.1149/1.2428700.
  * [30] J. B. Austin, Edgar collins bain, Biographical Memoirs 49 (1978) 25.
URL www.nasonline.org/publications/biographical-memoirs/memoir-pdfs/bain-
edgar-c.pdf

  * [31] H. Martin, P. Amoako-Yirenkyi, A. Pohjonen, N. K. Frempong, J. Komi, M. Somani, Statistical modeling for prediction of cct diagrams of steels involving interaction of alloying elements, Metallurgical and Materials Transactions B (2020) 1–13doi:10.1007/s11663-020-01991-w.
  * [32] H. Martin, A statistical and quantum mechanical model of combined alloying elements for phase transformation, Theses, Kwame Nkrumah University of Science and Technology, Kumasi, Ghana (2019).
  * [33] M. Henry, N. Oswald, B. C. D. Andrew, Improve the ductility of locally manufactured steel rods by tempering, in: 27th Biennial Conference Proceeding, Journal of the Ghana Science Association, 2011.
  * [34] N. Medvedeva, Y. N. Gornostyrev, A. Freeman, Solid solution softening and hardening in the group-v and group-vi bcc transition metals alloys: First principles calculations and atomistic modeling, Physical Review B 76 (21) (2007) 212104.  doi:10.1103/PhysRevB.76.212104.
  * [35] Y. Zhao, G. Lu, Qm/mm study of dislocation-hydrogen/helium interactions in $\alpha$-fe, Modelling and Simulation in Materials Science and Engineering 19 (6) (2011) 065004.  doi:10.1088/0965-0393/19/6/065004.
  * [36] M. Itakura, H. Kaburaki, M. Yamaguchi, T. Okita, The effect of hydrogen atoms on the screw dislocation mobility in bcc iron: a first-principles study, Acta Materialia 61 (18) (2013) 6857–6867.  doi:10.1016/j.actamat.2013.07.064.
  * [37] L. Romaner, V. Razumovskiy, R. Pippan, Core polarity of screw dislocations in fe–co alloys, Philosophical Magazine Letters 94 (6) (2014) 334–341.  doi:10.1080/09500839.2014.904055.
  * [38] L. Ventelon, B. Lüthi, E. Clouet, L. Proville, B. Legrand, D. Rodney, F. Willaime, Dislocation core reconstruction induced by carbon segregation in bcc iron, Physical Review B 91 (22) (2015) 220102.  doi:10.1103/PhysRevB.91.220102.
  * [39] Y.-J. Hu, M. R. Fellinger, B. G. Bulter, Y. Wang, K. A. Darling, L. J. Kecskes, D. R. Trinkle, Z.-K. Liu, Solute-induced solid-solution softening and hardening in bcc tungsten, Acta Materialia 141 (2017) 304–316.  doi:10.1016/j.actamat.2017.09.019.
  * [40] Y. Zhao, J. Marian, Direct prediction of the solute softening-to-hardening transition in w–re alloys using stochastic simulations of screw dislocation motion, Modelling and Simulation in Materials Science and Engineering 26 (4) (2018) 045002.  doi:10.1088/1361-651X/aaaecf.
  * [41] B. Lüthi, L. Ventelon, C. Elsässer, D. Rodney, F. Willaime, First principles investigation of carbon-screw dislocation interactions in body-centered cubic metals, Modelling and Simulation in Materials Science and Engineering 25 (8) (2017) 084001.  doi:10.1088/1361-651X/aa88eb.
  * [42] J. Svoboda, F. Fischer, Anisotropy of interstitial diffusion in bcc-crystals due to stress-induced unequal occupancies of different types of sites, International Journal of Solids and Structuresdoi:10.1016/j.ijsolstr.2018.05.023.
  * [43] Y. Zhao, L. Dezerald, J. Marian, Electronic structure calculations of oxygen atom transport energetics in the presence of screw dislocations in tungsten, Metals 9 (2) (2019) 252.  doi:10.3390/met9020252.
  * [44] B. Lüthi, L. Ventelon, D. Rodney, F. Willaime, Attractive interaction between interstitial solutes and screw dislocations in bcc iron from first principles, Computational Materials Science 148 (2018) 21–26.  doi:10.1016/j.commatsci.2018.02.016.
  * [45] Y. Wang, X. Wang, Q. Li, B. Xu, W. Liu, Atomistic simulations of carbon effect on kink-pair energetics of bcc iron screw dislocations, Journal of Materials Science 54 (15) (2019) 10728–10736.  doi:10.1007/s10853-019-03564-y.
  * [46] J. R. Stephens, Effects of interstitial impurities on the low-temperature tensile properties of tungsten, Tech. rep., National Aeronautics and Space Administration Cleveland OH Lewis Research Center (1964).
URL https://apps.dtic.mil/docs/citations/ADA396979

  * [47] M. Ulitchny, R. Gibala, The effects of interstitial solute additions on the mechanical properties of niobium and tantalum single crystals, Journal of the Less Common Metals 33 (1) (1973) 105–116.  doi:10.1016/0022-5088(73)90061-1.
  * [48] X. Jiang, M. Wang, K. Schmidt, E. Dunlop, J. Haupt, W. Gissler, Elastic constants and hardness of ion-beam-sputtered tin x films measured by brillouin scattering and depth-sensing indentation, Journal of applied physics 69 (5) (1991) 3053–3057.  doi:10.1063/1.348963.
  * [49] S.-H. Jhi, S. G. Louie, M. L. Cohen, J. Ihm, Vacancy hardening and softening in transition metal carbides and nitrides, Physical Review Letters 86 (15) (2001) 3348.  doi:10.1103/PhysRevLett.86.3348.
  * [50] X.-J. Chen, V. V. Struzhkin, Z. Wu, M. Somayazulu, J. Qian, S. Kung, A. N. Christensen, Y. Zhao, R. E. Cohen, H.-k. Mao, et al., Hard superconducting nitrides, Proceedings of the National Academy of Sciences 102 (9) (2005) 3198–3201.  doi:10.1073/pnas.0500174102.
  * [51] Z. Wu, X.-J. Chen, V. V. Struzhkin, R. E. Cohen, Trends in elasticity and electronic structure of transition-metal nitrides and carbides from first principles, Physical Review B 71 (21) (2005) 214103.  doi:10.1103/PhysRevB.71.214103.
  * [52] R. Hickman, T. McKechnie, Advanced materials and processes for boost phase nozzles, in: 40th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, 2004, p. 3385.  doi:10.2514/6.2004-3385.
  * [53] P. Thakre, V. Yang, Chemical erosion of refractory-metal nozzle inserts in solid-propellant rocket motors, Journal of Propulsion and Power 25 (1) (2009) 40–50.  doi:10.2514/1.37922.
  * [54] J.-Y. Yan, A. Ruban, Configurational thermodynamics of c in body-centered cubic/tetragonal fe: A combined computational study, Computational Materials Science 147 (2018) 293–303.  doi:10.1016/j.commatsci.2018.02.024.
  * [55] P. Hohenberg, W. Kohn, Inhomogeneous electron gas, Physical review 136 (3B) (1964) B864.  doi:10.1103/PhysRev.136.B864.
URL https://doi.org/10.1103/PhysRev.136.B864

  * [56] W. Kohn, L. J. Sham, Self-consistent equations including exchange and correlation effects, Physical review 140 (4A) (1965) A1133.
URL https://doi.org/10.1103/PhysRev.140.A1133

  * [57] G. Kresse, J. Furhmuller, Software vasp, vienna (1999); g. kresse, j. hafner, Phys. Rev. B 47 (R558).
  * [58] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Physical review B 54 (16) (1996) 11169.
URL https://doi.org/10.1103/PhysRevB.54.11169

  * [59] G. Kresse, J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational materials science 6 (1) (1996) 15–50.  doi:10.1016/0927-0256(96)00008-0.
  * [60] P. E. Blöchl, Projector augmented-wave method, Physical review B 50 (24) (1994) 17953.
URL https://doi.org/10.1103/PhysRevB.50.17953

  * [61] G. Kresse, D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Physical Review B 59 (3) (1999) 1758.
URL https://doi.org/10.1103/PhysRevB.59.1758

  * [62] A. Dewaele, M. Torrent, P. Loubeyre, M. Mezouar, Compression curves of transition metals in the mbar range: Experiments and projector augmented-wave calculations, Physical Review B 78 (10) (2008) 104102.
  * [63] B. J. Min, Study of the electronic and the structural properties of small molybdenum clusters via projector augmented wave pseudopotential calculations, Journal of the Korean Physical Society 66 (2) (2015) 209–213.
  * [64] L. Liang, Ab initio simulation of extended defects of &-Ti in presence of interstitial atoms H & O, Theses, Université Paris-Saclay (2016).
URL https://pastel.archives-ouvertes.fr/tel-01355132

  * [65] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Physical review letters 77 (18) (1996) 3865.  doi:10.1103/PhysRevLett.77.3865.
URL https://doi.org/10.1103/PhysRevLett.77.3865

  * [66] J. P. Perdew, M. Ernzerhof, A. Zupan, K. Burke, Nonlocality of the density functional for exchange and correlation: Physical origins and chemical consequences, The Journal of chemical physics 108 (4) (1998) 1522–1531.  doi:10.1063/1.475524.
  * [67] A. Zupan, K. Burke, M. Ernzerhof, J. P. Perdew, Distributions and averages of electron density parameters: Explaining the effects of gradient corrections, The Journal of chemical physics 106 (24) (1997) 10184–10193.  doi:10.1063/1.474101.
  * [68] A. Baldereschi, Mean-value point in the brillouin zone, Physical Review B 7 (12) (1973) 5212.
  * [69] H. J. Monkhorst, J. D. Pack, Special points for brillouin-zone integrations, Physical review B 13 (12) (1976) 5188.  doi:10.1103/PhysRevB.13.5188.
URL https://doi.org/10.1103/PhysRevB.13.5188

  * [70] M. Methfessel, A. Paxton, High-precision sampling for brillouin-zone integration in metals, Physical Review B 40 (6) (1989) 3616.
URL https://doi.org/10.1103/PhysRevB.40.3616

  * [71] D. D. Johnson, Modified broyden’s method for accelerating convergence in self-consistent calculations, Physical Review B 38 (18) (1988) 12807.
  * [72] W. P. Davey, Precision measurements of the lattice constants of twelve common metals, Physical Review 25 (6) (1925) 753.
  * [73] C. Kittel, P. McEuen, P. McEuen, Introduction to solid state physics, Vol. 8, Wiley New York, 1996.
  * [74] K. Hermann, Crystallography and Surface Structure: An Introduction for Surface Scientists and Nanoscientists, 2nd Edition, Wiley-VCH Verlag GmbH & Co. KGaA, 2017.
  * [75] E. I. Isaev, S. I. Simak, I. Abrikosov, R. Ahuja, Y. K. Vekilov, M. Katsnelson, A. Lichtenstein, B. Johansson, Phonon related properties of transition metals, their carbides, and nitrides: A first-principles study, Journal of applied physics 101 (12) (2007) 123519.  doi:10.1063/1.2747230.
  * [76] D. E. Gray, et al., American Institute of Physics handbook: Section editors: Bruce H. Billings [and others] Coordinating editor: Dwight E. Gray, McGraw-Hill Companies, 1972.
  * [77] M. Delheusy, X-ray investigation of nb/o interfaces, Theses, Universite Paris-Sud XI (2008).  doi:10.18419/opus-6994.
  * [78] T. Amriou, B. Bouhafs, H. Aourag, B. Khelifa, S. Bresson, C. Mathieu, Fp-lapw investigations of electronic structure and bonding mechanism of nbc and nbn compounds, Physica B: Condensed Matter 325 (2003) 46–56.  doi:10.1016/S0921-4526(02)01429-1.
  * [79] P. Haas, F. Tran, P. Blaha, Calculation of the lattice constant of solids with semilocal functionals, Physical Review B 79 (8) (2009) 085104.
  * [80] D. Connétable, First-principles study of transition metal carbides, Materials Research Express 3 (12) (2016) 126502.  doi:10.1088/2053-1591/3/12/126502.
  * [81] T. Shimanouchi, Tables of molecular vibrational frequencies part 5, Journal of Physical and Chemical Reference Data 1 (1) (1972) 189–216.
  * [82] J. Gutow, A. Herráez, Nitrogen gas (n2), Tech. rep., University of Wisconsin Oshkosh (UWOSH) (2009).
URL
http://www.uwosh.edu/faculty_staff/gutow/Chem_371_S09/OASI%20Web%20Pages/N2%20Calculations/N2%20Calculations/N2%20Calculations.html

  * [83] A. Weber, E. A. McGinnis, The raman spectrum of gaseous oxygen, Journal of Molecular Spectroscopy 4 (1-6) (1960) 195–200.  doi:10.1016/0022-2852(60)90081-3.
  * [84] W. McGrath, A. R. J. P. Ubbelohde, The influence of molecular flexibility on the transfer of vibrational energy in hydrocarbons, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 227 (1168) (1954) 1–9.  doi:10.1098/rspa.1954.0275.


# Magnetogenesis from early structure formation due to Yukawa forces.

Ruth Durrer  Department of Theoretical Physics, Université de Genève, Quai E.
Ansermet 24 , Genève, 1211, Switzerland.  Alexander Kusenko  Department of
Physics and Astronomy, University of California, Los Angeles
Los Angeles, California, 90095-1547, USA  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), UTIAS
The University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Theoretical Physics
Department, CERN, 1211 Geneva 23, Switzerland

###### Abstract

Yukawa interactions can mediate relatively long-range attractive forces
between fermions in the early universe. Such a globally attractive interaction
creates an instability that can result in the growth of structure in the
affected species even during the radiation dominated era. The formation and
collapse of fermionic microhalos can create hot fireballs at the sites of the
collapsing halos which inject energy into the cosmic plasma. We show that the
injected energy can be partially converted into primordial magnetic fields and
we estimate the correlation scale and the power spectrum of these fields. We
show that they may be the seeds of the observed astrophysical magnetic fields.

The origin of astrophysical magnetic fields is not well understood. The
observed fields in galaxies and clusters can be explained by the dynamo
amplification Ruzmaikin  _et al._ (1988) of seed fields produced in the early
universe or at later time during structure formation Durrer and Neronov
(2013). However, there are convincing arguments that small magnetic fields of
the order of $10^{-16}$ Gauss are also present in voids Neronov and Vovk
(2010); Essey  _et al._ (2011). Such fields are very difficult to produce at
late times and are therefore most certainly primordial. Future observations
will probe the structure of the magnetic fields in the intergalactic voids,
where they are expected to remain closely related to the primordial seed
fields. Combined with a theoretical understanding of magnetogenesis, this may
open a new window on the early universe.

In this Letter we identify a new class of scenarios in which magnetic fields
with astrophysically relevant correlation lengths are generated. It is based
on a recently discovered phenomenon of primordial structure formation in some
specific species interacting via Yukawa forces Amendola  _et al._ (2018);
Savastano  _et al._ (2019); Casas  _et al._ (2016); Flores and Kusenko
(2021a); Domènech and Sasaki (2021); Flores and Kusenko (2021b).

While significant gravitational growth of structure is possible during the
matter dominated era only, there is a possibility for structure formation in
some specific species interacting via other long-range attractive forces, even
in a radiation dominated universe Amendola  _et al._ (2018); Savastano  _et
al._ (2019); Casas  _et al._ (2016); Flores and Kusenko (2021a); Domènech and
Sasaki (2021); Flores and Kusenko (2021b). One simple example is the standard
model augmented by a heavy fermion $\psi$ and a light or massless boson
$\phi$, interacting via a Yukawa force. The attractive, relatively long-range
Yukawa forces create an instability, which is similar to the gravitational
instability, only much stronger. The growth and collapse of the halos of
fermions can lead to bound states Wise and Zhang (2014); Gresham  _et al._
(2018); Gresham and Zurek (2019), or, thanks to radiative cooling by the same
Yukawa interactions, it can lead to formation of primordial black holes Flores
and Kusenko (2021a, b). We shall call this possibility scenario 1.
Alternatively, the halos can disappear due to fermion annihilations in the
dense halos. We shall call this possibility scenario 2.

If the particles $\phi$ and $\psi$ interact with standard model particles, the
formation, collapse, and decay of the halos can locally inject energy into the
ambient plasma. In this Letter we show that such inhomogeneous energy
injection into the cosmic plasma creates suitable conditions for
magnetogenesis Durrer and Neronov (2013).

Let us illustrate the scenario on a simple model that involves one fermion and
one boson interacting via a Yukawa coupling:

$\mathcal{L}\supset\frac{1}{2}m_{\phi}^{2}{\phi}^{2}-y{\phi}\bar{\psi}\psi+\cdots.$
(1)

As discussed in Refs. Amendola  _et al._ (2018); Savastano  _et al._ (2019);
Casas  _et al._ (2016); Flores and Kusenko (2021a); Domènech and Sasaki
(2021); Flores and Kusenko (2021b); Flores  _et al._ (2022), the growth of
structures in $\psi$ particles due to the Yukawa instability, once it begins,
proceeds very rapidly and reaches the nonlinear regime within a Hubble time
$H^{-1}$. Each halo collapses and virializes, but scalar bremsstrahlung
($\psi\psi\rightarrow\psi\psi\phi$) provides radiative cooling and removes
energy from the halo, allowing it to collapse further. The outcome of this
process depends on the presence of fermion number asymmetry in the dark
sector. In the absence of an asymmetry ($n_{\psi}=n_{\bar{\psi}}$),
annihilations $\bar{\psi}\psi\rightarrow\phi\phi$ proceed at a growing rate
and eventually decimate the halo (scenario 2). Formation and disappearance of
halos in this case can set the right conditions for baryogenesis Flores  _et
al._ (2022). Alternatively, if the dark fermions develop an asymmetry, which
can be co-generated with the baryon asymmetry (as in the models of asymmetric
dark matter Petraki and Volkas (2013); Zurek (2014)), the outcome of halo
collapse is formation of primordial black holes (PBH), which can be dark
matter (scenario 1). The dark matter abundance is natural in the case of GeV
scale particles with an asymmetry comparable to the baryon asymmetry of the
universe:

$\Omega_{\rm PBH}=0.2\,{\Omega_{\rm
DM}}\frac{m_{\psi}}{m_{p}}\frac{\eta_{\psi}}{\eta_{\rm B}}={\Omega_{\rm
DM}}\left(\frac{m_{\psi}}{5\,{\rm
GeV}}\right)\left(\frac{\eta_{\psi}}{10^{-10}}\right).$ (2)

Even in the case of co-generation of the asymmetry in the baryonic and the
dark sector, the two parameters, $\eta_{B}$ and $\eta_{\psi}$ may differ by a
few orders of magnitude due to the differences in the temperatures and the
numbers of degrees of freedom in the visible and the dark sectors,
respectively Petraki and Volkas (2013); Zurek (2014). The formation of
structures in $\psi$ particles is possible after the fermions kinetically
decouple from the plasma, which happens at temperature $T_{*}\sim
10^{-2}m_{\psi}$ Graesser  _et al._ (2011); Flores and Kusenko (2021a). In
what follows, we will assume the following representative values of the
parameters, which are consistent with dark matter in the form of PBHs:

$\displaystyle T_{*}$ $\displaystyle\sim$ $\displaystyle 1-100\ {\rm MeV},$
(3) $\displaystyle m_{\psi}$ $\displaystyle\sim$ $\displaystyle 0.1-10\ {\rm
GeV}.$ (4)

It is reasonable to assume that the dark sector $\\{\psi,\phi\\}$ is weakly
coupled to the standard model via some higher-dimensional operators or small
kinetic mixing, which is, in fact, necessary for the co-generation of the
fermion asymmetry Petraki and Volkas (2013); Zurek (2014). This implies that
during the halo collapse and emission of scalar particles, a small fraction of
the halo energy is transmitted to the cosmic plasma. In scenario 1, where the
halos finally form the dark matter, the energy transmitted to the plamsa is a
fraction of the dark matter density, $\epsilon_{1}\rho_{DM}$, while in
scenario 2 it is a faction of the radiation density, $\epsilon_{2}\rho_{r}$.
This energy injection results in creation of a fireball at the site of each
halo. As the fireballs expand, since the Reynolds numbers of the plasma are
very high, a significant fraction of this inhomogeneously deposited energy
will be converted into turbulence. We denote this fraction $f$. We thus
parameterize the energy in turbulence as a fraction of the dark matter energy:
$\rho_{\rm turb}=\epsilon_{1}f\,\rho_{\rm DM}=\epsilon\rho_{\rm DM}$ in
scenario1, while $\rho_{\rm turb}=\epsilon_{2}f\,\rho_{r}=\epsilon\rho_{r}$ in
scenario 2.

As electrons and positrons are still relativistic at $T_{*}$, when turbulence
is generated, we can assume the turbulence to develop within a Hubble time
into Kolmogorov turbulence with a correlation scale of the order of the Hubble
scale.

During the radiation era,

$\displaystyle H(T_{*})$ $\displaystyle=$
$\displaystyle\frac{1}{2\tau_{*}}~{}\simeq~{}5\sqrt{g_{\rm
eff\,*}}\left(\frac{T_{*}}{1{\rm MeV}}\right)^{2}{\rm sec}^{-1}$ (5)
$\displaystyle\lambda_{*}$ $\displaystyle=$
$\displaystyle\eta\frac{1+z_{*}}{H(T_{*})}~{}\simeq~{}\eta\frac{1.5}{\sqrt{g_{\rm
eff\,*}}}10^{11}\left(\frac{1{\rm MeV}}{T_{*}}\right){\rm sec}$ (6)
$\displaystyle\simeq$ $\displaystyle\frac{1.7\eta}{\sqrt{g_{\rm
eff\,*}}}\left(\frac{1{\rm MeV}}{T_{*}}\right){\rm kpc}\,.$

Here $\tau$ is cosmic time and $\lambda_{*}$ is the comoving correlation scale
at the time $\tau_{*}$, at temperature $T_{*}$ and $g_{\rm eff\,*}$ is the
number of relativistic degrees of freedom at $T_{*}$. The speed of light is
set to $c=1$ and we normalize the scale factor to one today, $a_{0}=1$. The
factor $\eta\lesssim 1$ is ratio $\lambda_{*}/H_{*}^{-1}$. We assume that
rapidly (i.e. within less than one Hubble time), on wave numbers $k$ with
$2\pi/\lambda_{*}=k_{*}\leq k\leq k_{d}$ a Kolmogorov spectrum is established.
The scale $k_{d}$ is a damping scale beyond which the turbulent motion is
converted into heat.

Furthermore, the relativistic plasma is compressible, $\partial
P/\partial\rho=1/3$, hence compressible turbulence with a white noise spectrum
on large scales, $k<k_{*}=2\pi/\lambda_{*}$ develops. The turbulent velocity
power spectrum of causally generated compressible turbulence can be
approximated as follows

$\displaystyle\langle v_{i}({\mathbf{k}})v_{j}({\mathbf{k}}^{\prime})\rangle$
$\displaystyle=$
$\displaystyle(2\pi)^{3}\delta({\mathbf{k}}-{\mathbf{k}}^{\prime})\delta_{ij}P_{K}(k)$
(7)

with

$\displaystyle P_{K}(k)$ $\displaystyle=$ $\displaystyle
v_{*}^{2}k_{*}^{-3}\times\left\\{\begin{array}[]{ll}1\,,&k<k_{*}\\\
(k/k_{*})^{-11/3}\,,&k_{*}\leq k<k_{d}\\\ 0\,,&k_{d}\leq
k\,.\end{array}\right.$ (11)

Even though Kolmogorov has derived the spectral index $-11/3$ only for non-
relativistic turbulence, it has been shown in numerical simulations that also
relativistic turbulence develops a (nearly) Kolmogorov spectrum Zhdankin  _et
al._ (2017); Comisso and Sironi (2018). The amplitude $v_{*}^{2}$ is
proportional to $\epsilon$.

The energy density in the turbulence is given by the integral of the velocity
power spectrum,

$\displaystyle\rho_{K}$ $\displaystyle\simeq$
$\displaystyle\frac{\rho_{r}}{2}\langle v(x)^{2}\rangle$ (12) $\displaystyle=$
$\displaystyle\frac{3\rho_{r}}{4\pi^{2}}\int
P_{K}(k)k^{2}dk=\frac{11}{8\pi^{2}}\rho_{r}v_{*}^{2}\,.$

Here $\rho_{r}$ is the energy in relativistic particles which dominates the
total energy density in the radiation era.

As the cosmic plasma is charged and has very high conductivity, the turbulence
becomes MHD turbulence and will also lead to the generation of magnetic
fields. The magnetic field spectrum is of the form

$\langle
B_{i}({\mathbf{k}})B_{j}^{*}({\mathbf{k}}^{\prime})\rangle=(2\pi)^{3}\delta({\mathbf{k}}-{\mathbf{k}}^{\prime})\left(\delta_{ij}-\hat{k}_{i}\hat{k}_{j}\right)P_{B}(k)$
(13)

with

$P_{B}(k)=B_{*}^{2}k_{*}^{-3}\times\left\\{\begin{array}[]{ll}(k/k_{*})^{2}\,,&k<k_{*}\\\
(k/k_{*})^{-11/3}\,,&k_{*}\leq k<k_{d}\\\ 0\,,&k_{d}\leq
k\,.\end{array}\right.$ (14)

The $k^{2}$ behavior on large scales is due to the fact that the spectrum must
be analytic for small $k$ as the correlation function in real space has
compact support Durrer and Caprini (2003) (there are no correlations on super-
horizon scales). Note that we always calculate the comoving magnetic field.
The true magnetic field is $B/a^{2}$. The energy density in this magnetic
field is given by $\rho_{B}$, where (we use Heavyside units so that $B$ is
canonically normalized)

$\displaystyle a^{4}\rho_{B}(\tau_{*})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\langle
B({\mathbf{x}})^{2}\rangle=\frac{1}{2\pi^{2}}\int_{0}^{\infty}P_{B}(k)k^{2}dk$
(15) $\displaystyle=$
$\displaystyle\frac{1/5+3/2}{2\pi^{2}}B_{*}^{2}=\frac{17}{20\pi^{2}}B_{*}^{2}\,.$

Typically we expect equi-partition, $\rho_{B}\simeq\rho_{K}$.

We want to express $\epsilon$ in terms of the magnetic field $B_{*}$ for both
scenarios. The ratio of the energy density in magnetic fields and in thermal
radiation of $g_{\rm eff}$ degrees of freedom at temperature $T_{*}$ is given
by

$\frac{\rho_{B}}{\rho_{r}}=\frac{a^{4}\rho_{B}}{a^{4}\rho_{r}}=\frac{\langle
B^{2}\rangle}{g_{\rm eff}\left(2\times 10^{-6}{\rm Gauss}\right)^{2}}\,.$ (16)

Therefore, for scenario 2

$\epsilon\simeq\frac{\rho_{B}}{\rho_{r}}\simeq
10^{-18}\left(\frac{B_{*}}{10^{-15}{\rm Gauss}}\right)^{2}\,.$ (17)

For scenario 1 we need the ratio of $\rho_{B}$ and $\rho_{DM}$. At the time
$\tau_{*}$, the Universe is radiation dominated, hence the fraction of the
dark matter density dumped into magnetic fields can be expressed as

$\displaystyle\frac{\rho_{B}(\tau_{*}}{\rho_{\rm DM}(\tau_{*})}$
$\displaystyle\simeq$
$\displaystyle\frac{\rho_{B}}{\rho_{r}}\frac{1+z_{*}}{1+z_{eq}}\simeq\frac{\rho_{B}}{\rho_{r}}\frac{T_{*}/T_{0}}{3200}$
(18) $\displaystyle\simeq$ $\displaystyle 1.3\times
10^{6}\left(\frac{T_{*}}{1{\rm MeV}}\right)\frac{\langle B^{2}\rangle}{g_{\rm
eff}\left(2\times 10^{-6}{\rm Gauss}\right)^{2}}$ $\displaystyle\simeq$
$\displaystyle 10^{-12}\left(\frac{T_{*}}{1{\rm
MeV}}\right)\left(\frac{B_{*}}{10^{-15}{\rm Gauss}}\right)^{2}\,,$

were $T_{0}=2.73$K is the present temperature of the Universe. Hence for
$B_{*}=10^{-15}$ Gauss one only needs $\epsilon\sim 10^{-18}$ in scenario 2.
In scenario 1, the value if $\epsilon$ also depends on the temperature. At
$T_{*}\sim 1$MeV one needs $\epsilon\sim 10^{-12}$, while the same amplitude
magnetic fields at $T_{*}\sim 10^{5}$GeV requires $\epsilon\sim 10^{-4}$.

If the magnetic fields would simply scale like $a^{-2}$, $B_{*}$ would be the
physical magnetic field today. However, MHD turbulence once generated is
freely decaying with an inverse cascade due to the white noise velocity power
spectrum. This leads to significant changes in both, the amplitude of the
magnetic field and its correlation scale.which can grow up to nearly 1Mpc.

At $T_{*}$ the situation is as discussed above. Let us study its evolution.

Equi-partition requests that $v_{*}\simeq v_{A}=\sqrt{\rho_{B}/\rho_{r}}$,
where $v_{A}$ is the Alfvén velocity. The comoving power spectra and the
correlations scale $k_{*}$ evolve with time (see Durrer and Neronov (2013)
where this is derived in detail). More precisely, $\lambda_{*}$ behaves as

$\lambda_{*}(t)\propto t^{\frac{2}{5}}\,.$ (19)

Here we use the fact that we consider compressible turbulence, where the
velocity has a white noise spectrum. Therefore, correlation scale grows
significantly as equi-partition requires the peak of the magnetic field power
spectrum to move to the left (see Fig. 1), provoking an inverse cascade. The
energy density evolves as

$a^{4}\rho_{B}\simeq a^{4}\rho_{K}\propto\frac{\overline{v}^{2}}{2}\propto
t^{-6/5}$ (20)

Here and above, $t$ is conformal time. Also this result is derived in detail
in Durrer and Neronov (2013) and has also been confirmed by numerical
simulations Hindmarsh  _et al._ (2014). Schematically the evolution looks as
shown in Fig. 1.

Figure 1: The evolution of the magnetic field spectrum in an MHD plasma with
compressible turbulence. The spectra evolve towards smaller $k$. Dashed lines
show the kinetic energy spectrum. In this schematic diagram the decay of the
peak amplitude is not taken into account (Figure adapted from Durrer and
Neronov (2013)).

This evolution continues until about recombination after which there is no
longer a charged plasma and the evolution (roughly) stops.

At late times we find a correlation scale of

$\displaystyle\lambda_{*}(t_{0})$ $\displaystyle\simeq$
$\displaystyle\lambda_{*}(t_{\rm rec})=\lambda_{*}(t_{*})\left(\frac{t_{\rm
rec}}{t_{*}}\right)^{2/5}$ (21) $\displaystyle\simeq$
$\displaystyle\lambda_{*}(t_{*})\left(\frac{z_{*}+1}{z_{\rm
rec}+1}\right)^{2/5}$ $\displaystyle\simeq$
$\displaystyle\lambda_{*}(t_{*})\left(\frac{(T_{*}/1{\rm
MeV})10^{10}}{2.35\times 1000}\right)^{2/5}$ $\displaystyle\simeq$
$\displaystyle 450\times\lambda_{*}(t_{*})\left(\frac{T_{*}}{1{\rm
MeV}}\right)^{2/5}$ $\displaystyle\simeq$ $\displaystyle 0.76{\rm
Mpc}\frac{\eta}{\sqrt{g_{\rm eff}}}\left(\frac{1{\rm
MeV}}{T_{*}}\right)^{3/5}\,.$

Correspondingly one finds for the evolved magnetic field

$\displaystyle B_{0}$ $\displaystyle\equiv$ $\displaystyle B_{*}(t_{0})\simeq
B_{*}(t_{\rm rec})\simeq B_{*}(t_{*})\left(\frac{t_{\rm
rec}}{t_{*}}\right)^{-3/5}$ $\displaystyle\simeq$ $\displaystyle
B_{*}(t_{*})\left(\frac{z_{*}+1}{z_{\rm dec}+1}\right)^{-3/5}$
$\displaystyle\simeq$ $\displaystyle B_{*}(t_{*})\left(\frac{(T_{*}/1{\rm
MeV})10^{10}}{2.35\times 1000}\right)^{-3/5}$ (22) $\displaystyle\simeq$
$\displaystyle 10^{-4}\times B_{*}(t_{*})\left(\frac{T_{*}}{1{\rm
MeV}}\right)^{-3/5}$

In this case, in order to have magnetic fields of about $10^{-15}$Gauss today,
we need $B_{*}(t_{*})\simeq 10^{-11}$ Gauss. This means that at time
$\tau_{*}$ roughly we need in general

Scenario 1 : $\displaystyle\epsilon=\frac{\rho_{B}(\tau_{*})}{\rho_{\rm
DM}(\tau_{*})}$ $\displaystyle\simeq$ $\displaystyle
10^{-6}\left(\frac{T_{*}}{1{\rm
MeV}}\right)^{11/5}\left(\frac{B_{0}}{10^{-16}{\rm Gauss}}\right)^{2}\,,$
Scenario 2 :
$\displaystyle\epsilon=\frac{\rho_{B}(\tau_{*})}{\rho_{r}(\tau_{*})}$
$\displaystyle\simeq$ $\displaystyle 10^{-12}\left(\frac{T_{*}}{1{\rm
MeV}}\right)^{6/5}\left(\frac{B_{0}}{10^{-16}{\rm Gauss}}\right)^{2}\,.$ (24)

Here $B_{0}$ is the (physical) magnetic field amplitude today.

In Fig. 2 we plot the required initial ratio of $\epsilon=\rho_{B}/\rho_{\rm
DM}$ for scenario 1 (top panel) and $\epsilon=\rho_{B}/\rho_{r}$ for scenario
2 (bottom panel) needed to obtain a field amplitude $B_{0}$ today for a given
explosion temperature $T_{*}$. For scenario 1, the two brightest colors with
$\epsilon\geq 1$ are clearly excluded. We must require $\epsilon<0.1$ A
present amplitude of about $10^{-14}$Gauss with a formation temperature of
$1$MeV would require $\epsilon=0.01$ which seems reasonable. Also an amplitude
of $10^{-16}$ to $10^{-15}$ Gauss at $T_{*}=10$MeV to $T_{*}=1$MeV can be
obtained with a modest value of $\epsilon\simeq 10^{-4}$. The results for
scenario 2 are much more optimistic. For all the injection temperatures
considered here, $0.1$MeV$\leq T_{*}\leq 100$MeV present magnetic fields of
nearly up to $10^{-12}$Gauss can be generated with $\epsilon<10^{-2}$. For a
present magnetic field of $10^{-16}$Gauss even a value of $\epsilon\leq
10^{-9}$ is sufficient.

Figure 2: The top panel shows the required ratio,
$\log\epsilon=\log\left(\frac{\rho_{B}(\tau_{*})}{\rho_{\rm
DM}(\tau_{*})}\right)$ as a function of temperature $T_{*}$ and final magnetic
field $B_{0}$ (the log’s are to base 10) for scenario 1.
The bottom panel shows the required ratio,
$\log\epsilon=\log\left(\frac{\rho_{B}(\tau_{*})}{\rho_{r}(\tau_{*})}\right)$
as a function of temperature $T_{*}$ and final magnetic field $B_{0}$ (the
log’s are to base 10) for scenario 2.

In this work we have found that explosions such as they are expected from
early structure formation in the radiation dominated Universe at $100$MeV
$>T_{*}>0.1$MeV can lead to magnetic fields with an amplitude of about
$10^{-16}$Gauss or more. The correlation scale of these fields can be of order
Mpc for $T_{*}=1$MeV and somewhat smaller/larger at higher/lower temperature.
A significantly lower temperature is most probably in conflict with other
observations such as primordial nucleosynthesis. In the first scenario, where
primordial black holes are formed which are the dark matter, explosions at a
much higher temperature generate magnetic fields that have too low amplitude
today due to dilution and turbulent decay. However, in the second scenario
where the dark sector fermions annihilate into radiation, magnetic fields with
amplitude up to $10^{-16}$Gauss can be generated even at $T_{*}\sim
10^{6}$GeV, albeit with much shorter correlation length. Note also, that the
important ingredient here was just that we have fast energy injection into the
cosmic plasma at some high temperature $T_{*}$, correlated over a scale close
to the Hubble scale $H^{-1}(T_{*})$, and that we can characterize the injected
energy as $\epsilon\rho_{DM}$ (scenario 1) or $\epsilon\rho_{r}$ (scenario 2).
The details of the energy injection mechanism are irrelevant for these
results.

If, indeed, the seed magnetic fields are produced in the early universe and
are amplified in galaxies and clusters, the fields deep in the voids are
expected to retain information about the original primordial magnetic power
spectrum. The upcoming gamma-ray observations with the Cherenkov Telescope
Array (CTA) can provide new information because the gamma-ray cascades from
blazars are sensitive to both the amplitude and the correlation lengths of the
magnetic fields in the intergalactic voids Neronov and Semikoz (2007); Neronov
and Vovk (2010); Essey  _et al._ (2011); Dermer  _et al._ (2011); Razzaque
_et al._ (2012). Thus future observations could open a new window on the early
universe and elucidate the processes that generated primordial magnetic seed
fields.

###### Acknowledgements.

We thank Guillem Domènech, Marcos Flores, Derek Inman, Andrii Neronov, and
Misao Sasaki for helpful discussions. R.D. acknowledges support from the Swiss
National Science Foundation, grant No. 200020 182044\. A.K. was supported by
the U.S. Department of Energy (DOE) Grant No. DE-SC0009937, by the Simons
Foundation Fellowship, by the World Premier International Research Center
Initiative (WPI), MEXT, Japan, by Japan Society for the Promotion of Science
(JSPS) KAKENHI grant No. JP20H05853, and by the UC Southern California Hub
with funding from the UC National Laboratories division of the University of
California Office of the President.

## References

  * Ruzmaikin  _et al._ (1988) A. A. Ruzmaikin, D. D. Sokolov, and A. M. Shukurov,  _Magnetic Fields of Galaxies_, Vol. 133 (Springer Verlag, Dortrecht, 1988).
  * Durrer and Neronov (2013) R. Durrer and A. Neronov, Astron. Astrophys. Rev. 21, 62 (2013), arXiv:1303.7121 [astro-ph.CO] .
  * Neronov and Vovk (2010) A. Neronov and I. Vovk, Science 328, 73 (2010), arXiv:1006.3504 [astro-ph.HE] .
  * Essey  _et al._ (2011) W. Essey, S. Ando, and A. Kusenko, Astropart. Phys. 35, 135 (2011), arXiv:1012.5313 [astro-ph.HE] .
  * Amendola  _et al._ (2018) L. Amendola, J. Rubio, and C. Wetterich, Phys. Rev. D 97, 081302 (2018), arXiv:1711.09915 [astro-ph.CO] .
  * Savastano  _et al._ (2019) S. Savastano, L. Amendola, J. Rubio, and C. Wetterich, Phys. Rev. D 100, 083518 (2019), arXiv:1906.05300 [astro-ph.CO] .
  * Casas  _et al._ (2016) S. Casas, V. Pettorino, and C. Wetterich, Phys. Rev. D 94, 103518 (2016), arXiv:1608.02358 [astro-ph.CO] .
  * Flores and Kusenko (2021a) M. M. Flores and A. Kusenko, Phys. Rev. Lett. 126, 041101 (2021a), arXiv:2008.12456 [astro-ph.CO] .
  * Domènech and Sasaki (2021) G. Domènech and M. Sasaki, JCAP 06, 030 (2021), arXiv:2104.05271 [hep-th] .
  * Flores and Kusenko (2021b) M. M. Flores and A. Kusenko, (2021b), arXiv:2108.08416 [hep-ph] .
  * Wise and Zhang (2014) M. B. Wise and Y. Zhang, Phys. Rev. D 90, 055030 (2014), [Erratum: Phys.Rev.D 91, 039907 (2015)], arXiv:1407.4121 [hep-ph] .
  * Gresham  _et al._ (2018) M. I. Gresham, H. K. Lou, and K. M. Zurek, Phys. Rev. D 97, 036003 (2018), arXiv:1707.02316 [hep-ph] .
  * Gresham and Zurek (2019) M. I. Gresham and K. M. Zurek, Phys. Rev. D 99, 083008 (2019), arXiv:1809.08254 [astro-ph.CO] .
  * Flores  _et al._ (2022) M. M. Flores, A. Kusenko, L. Pearce, and G. White, (2022), arXiv:2208.09789 [hep-ph] .
  * Petraki and Volkas (2013) K. Petraki and R. R. Volkas, Int. J. Mod. Phys. A28, 1330028 (2013), arXiv:1305.4939 [hep-ph] .
  * Zurek (2014) K. M. Zurek, Phys. Rept. 537, 91 (2014), arXiv:1308.0338 [hep-ph] .
  * Graesser  _et al._ (2011) M. L. Graesser, I. M. Shoemaker, and L. Vecchi, JHEP 10, 110 (2011), arXiv:1103.2771 [hep-ph] .
  * Zhdankin  _et al._ (2017) V. Zhdankin, G. R. Werner, D. A. Uzdensky, and M. C. Begelman, Phys. Rev. Lett. 118, 055103 (2017), arXiv:1609.04851 [physics.plasm-ph] .
  * Comisso and Sironi (2018) L. Comisso and L. Sironi, Phys. Rev. Lett. 121, 255101 (2018), arXiv:1809.01168 [astro-ph.HE] .
  * Durrer and Caprini (2003) R. Durrer and C. Caprini, JCAP 11, 010 (2003), arXiv:astro-ph/0305059 .
  * Hindmarsh  _et al._ (2014) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Phys. Rev. Lett. 112, 041301 (2014), arXiv:1304.2433 [hep-ph] .
  * Neronov and Semikoz (2007) A. Neronov and D. V. Semikoz, JETP Lett. 85, 473 (2007), arXiv:astro-ph/0604607 .
  * Dermer  _et al._ (2011) C. D. Dermer, M. Cavadini, S. Razzaque, J. D. Finke, J. Chiang, and B. Lott, Astrophys. J. Lett. 733, L21 (2011), arXiv:1011.6660 [astro-ph.HE] .
  * Razzaque  _et al._ (2012) S. Razzaque, C. D. Dermer, and J. D. Finke, Astrophys. J. 745, 196 (2012), arXiv:1110.0853 [astro-ph.HE] .