# How does dark matter stabilize disc galaxies?

K. Aditya
Indian Institute of Astrophysics, Koramangala, Bengaluru 560 034, INDIA
E-mail<EMAIL_ADDRESS>
###### Abstract

The study presents a theoretical framework for understanding the role of dark
matter on the stability of the galactic disc. We model the galaxy as a two-
component system consisting of stars and gas in equilibrium with an external
dark matter halo. We derive the equations governing the growth of
perturbations and obtain a stability criterion that connects the potential of
the dark matter halo and the gas fraction with the stability levels of the
galaxy. We find that a two-component disc is more susceptible to the growth of
gravitational instabilities than individual components, particularly as gas
fractions increase. However, the external field, due to the dark matter halo,
acts as a stabilizing agent and increases the net stability levels even in the
presence of a cold gas component. We apply the stability criterion to models
of the Milky Way, low surface brightness galaxies, and baryon-dominated cold
rotating disc galaxies observed in the early universe. Our results show that
the potential due to the dark matter halo plays a significant role in
stabilizing nearby galaxies, such as the Milky Way, and low surface brightness
galaxies, which would otherwise be prone to local gravitational instabilities.
However, we find that the baryon-dominated cold disc galaxies observed in the
early universe remain susceptible to the growth of local gravitational
instabilities despite the stabilizing effect of the dark matter halo.

###### keywords:

hydrodynamics-instabilities, galaxies:kinematics and dynamics,
galaxies:structure, galaxies:star formation, Galaxy:evolution

##  1 Introduction

Gravitational instabilities are fundamental processes that drive the evolution
of the galaxy. It provides important clues for understanding how gas in the
galaxies is converted into stars (Wang & Silk, 1994; Dopita & Ryder, 1994;
Pandey & Van De Bruck, 1999; Krumholz & Burkert, 2010; Forbes et al., 2014),
and how non-axisymmetric structures like bars and spiral arms form in galaxies
(Goldreich & Lynden-Bell, 1965; Toomre, 1977; Iye, 1978; Kalnajs, 1983; Lin &
Shu, 1987; Sellwood, 2000, 2012). One of the simplest diagnostics for
accessing the stability of the galactic disc against the growth of
axisymmetric gravitational instabilities was proposed by Toomre (1964). It
measures the competing effect of self-gravity, which tries to destabilize the
disc, and the stabilizing effect of the differential rotation and the random
velocity dispersion. The balance between the stabilizing agents, i.e.,
differential rotation and random velocity dispersion, and the destabilizing
agent, i.e., the self-gravity, is classically quantified by the stability
criterion proposed by Toomre (1964):

$q=\frac{\kappa\sigma}{\pi G\Sigma}.$ (1)

In the above equation, $\kappa$ is the epicyclic frequency, $\Sigma$ is the
mass surface density, and $\sigma$ is the radial velocity dispersion, where
$q>1$ is the condition for stability of the disc against axisymmetric
perturbations. The stability criterion proposed by Toomre (1964) has been
modified to include the self-gravity of both stars and gas by Jog & Solomon
(1984a); Wang & Silk (1994); Elmegreen (1995); Jog (1996); Rafikov (2001);
Romeo & Wiegert (2011) and finally Romeo & Falstad (2013) derive a N-component
stability parameter to quantify the stability of gravitationally coupled
multiple stellar and gaseous discs. The stability parameter has also been
modified further to include the physical processes like the effects of the
turbulence (Hoffmann & Romeo, 2012; Agertz et al., 2015) and the three-
dimensional structure of ISM (Meidt, 2022).
The two-component model for studying the stability of disc galaxy against the
growth of local axisymmetrical gravitational instabilities was envisaged by
Jog & Solomon (1984a). In the two-component model, stars and gas in the
galactic disk are modeled as two isothermal fluids that interact
gravitationally with each other. One of the components in the two-component
model resembles the interstellar medium (ISM) with smaller values of the
velocity dispersion, and the other resembles the stellar component with higher
velocity dispersion. The approach has been used extensively to study the role
of the cold ISM in driving the instabilities in galactic disc (Jog & Solomon,
1984b) and for studying the stability of the Galactic disc by (Jog, 1996).
Further, Rafikov (2001) presents the stability criterion for a disc consisting
of multiple isothermal components. Each component is categorized as either
collisional, such as the ISM, or collisionless, such as stars. The results
obtained by Rafikov (2001) for the stability using a collisionless treatment
for stars and collisional approach for the ISM are comparable to the results
obtained Jog & Solomon (1984a) and Elmegreen (1995). The two-component
stability parameter is a valuable diagnostic for understanding if the
stability levels are driven by stars or by gas; for example, see Romeo &
Mogotsi (2017) and Romeo & Fathi (2016).
The stability criterion in the literature considers the self-gravity of the
gas and stars. However, it does not consider the role of dark matter in
driving the gravitational instabilities. The initial effort to incorporate the
influence of a dark matter halo on the stability of a single-component disc
was undertaken by Jog (2014). In this work, we present the conditions for
appraising the stability of the gravitationally coupled two-component disc
consisting of stars and gas in equilibrium with an external dark matter halo.
The differential equations governing the growth rate of perturbations are
derived by considering the two-component disc in equilibrium with the external
force field of the dark matter halo. Each component is specified by its
velocity dispersion, surface density, and angular frequency, but the system is
under the influence of the force field of the dark matter halo. We show that
the governing equations for the growth of instabilities resemble a wave
equation with extra terms. We use plane wave ansatz to derive the dispersion
relation for the gravitationally coupled two-component system in equilibrium
with the dark matter halo and obtain a simple stability criterion. The
stability criterion presented in this work explicitly quantifies the
contribution of dark matter to the overall stability levels. It can be used to
explore the role of dark matter in regulating various physical processes
within the galactic disc where gravitational instabilities are important.
The paper is organized as follows: in §2, we will formulate the basic
equations and derive the governing differential equations. We will derive the
dispersion relation and stability criterion in §3 and §4. We finally present
the results in §5 and discuss the applications of the stability criterion in
§6 $\&$ §7, and conclude in §8.

##  2 Formulation and derivation of basic equations

We consider a coaxial and coplanar thin disc comprising stars and gas, which
interact with each other gravitationally. The two-component disc is supported
by random pressure and rotation, and the system is in equilibrium with a
constant external force field of the dark matter halo. The problem is
described in the galactic cylindrical coordinate system $(R,\theta,z)$. We
start with the basic hydrodynamic equations in which the external force field
of the dark matter halo is in equilibrium with the two-component disc. We then
introduce small perturbations in the basic equations and derive the dynamic
equations governing the evolution of the perturbed quantities. The Force
equation, continuity equation, and the Poissons equation for a thin disc in
equilibrium with an external potential $\Phi_{ext}$ are:

$\Sigma_{i}\frac{\partial\boldsymbol{V_{i}}}{\partial
t}+\Sigma_{i}(\boldsymbol{V_{i}}.\boldsymbol{\nabla})\boldsymbol{V_{i}}=-\boldsymbol{\nabla}P_{i}-\Sigma_{i}\boldsymbol{\nabla}(\Phi_{s}+\Phi_{g})-\Sigma_{i}\boldsymbol{\nabla}\Phi_{ext},$
(2)

$\frac{\partial\Sigma_{i}}{\partial
t}+\boldsymbol{\nabla}.(\Sigma_{i}\boldsymbol{V_{i}})=0,$ (3)

$\nabla^{2}(\Phi_{s}+\Phi_{g})=4\pi G(\Sigma_{s}+\Sigma_{g})\delta(z).$ (4)

The above equations, when expressed in cylindrical coordinates, supplemented
with an isothermal equation of state $P_{i}=\Sigma_{i}c^{2}_{i}$ read:

$\Sigma_{i}\frac{\partial u_{i}}{\partial t}+\Sigma_{i}u_{i}\frac{\partial
u_{i}}{\partial R}+\frac{v_{i}\Sigma_{i}}{R}\frac{\partial
u_{i}}{\partial\theta}-\frac{\Sigma_{i}v^{2}_{i}}{R}=-c_{i}^{2}\frac{\partial\Sigma_{i}}{\partial
R}-\Sigma_{i}\frac{\partial(\Phi_{s}+\Phi_{g})}{\partial
R}-\Sigma_{i}\frac{\partial\Phi_{ext}}{\partial R},$ (5)

$\Sigma_{i}\frac{\partial v_{i}}{\partial t}+\Sigma_{i}u_{i}\frac{\partial
v_{i}}{\partial R}+\frac{\Sigma_{i}v_{i}}{R}\frac{\partial
v_{i}}{\partial\theta}+\frac{\Sigma_{i}v_{i}u_{i}}{R}=-\frac{c_{i}^{2}}{R}\frac{\partial\Sigma_{i}}{\partial\theta}-\frac{\Sigma_{i}}{R}\frac{\partial(\Phi_{s}+\Phi_{g})}{\partial\theta},$
(6)

$\frac{\partial\Sigma_{i}}{\partial
t}+\frac{1}{R}\frac{\partial(Ru_{i}\Sigma_{i})}{\partial
R}+\frac{v_{i}}{R}\frac{\partial\Sigma_{i}}{\partial\theta}+\frac{\Sigma_{i}}{R}\frac{\partial
v_{i}}{\partial\theta}=0,$ (7)

$\frac{1}{R}\frac{\partial}{\partial
R}\left(R\frac{\partial(\Phi_{s}+\Phi_{g})}{\partial
R}\right)+\frac{\partial^{2}(\Phi_{s}+\Phi_{g})}{\partial
z^{2}}+\frac{1}{R^{2}}\frac{\partial^{2}(\Phi_{s}+\Phi_{g})}{\partial\theta^{2}}=4\pi
G(\Sigma_{s}+\Sigma_{g})\delta(z).$ (8)

In the above equations, ${}^{\prime}i^{\prime}$ is used to index
${}^{\prime}stars^{\prime}$ and ${}^{\prime}gas^{\prime}$, $u_{i}$ and $v_{i}$
are the velocity components in the radial and the tangential directions
respectively, $\Sigma_{i}$ and $\Phi_{i}$ are the surface density and the
gravitational potential associated with the stellar and the gas disc,
respectively, and $c_{i}$ is the velocity dispersion of each component.
Assuming the disc is axisymmetric, the above equations can be written as

$\Sigma_{i}\frac{\partial u_{i}}{\partial t}+\Sigma_{i}u_{i}\frac{\partial
u_{i}}{\partial
R}-\frac{\Sigma_{i}v^{2}_{i}}{R}=-c_{i}^{2}\frac{\partial\Sigma_{i}}{\partial
R}-\Sigma_{i}\frac{\partial(\Phi_{s}+\Phi_{g})}{\partial
R}-\Sigma_{i}\frac{\partial\Phi_{ext}}{\partial R},$ (9)

$\Sigma_{i}\frac{\partial v_{i}}{\partial t}+\Sigma_{i}u_{i}\frac{\partial
v_{i}}{\partial R}+\frac{\Sigma_{i}v_{i}u_{i}}{R}=0,$ (10)

$\frac{\partial\Sigma_{i}}{\partial
t}+\frac{1}{R}\frac{\partial(Ru_{i}\Sigma_{i})}{\partial R}=0,$ (11)

$\frac{1}{R}\frac{\partial}{\partial
R}\left(R\frac{\partial(\Phi_{s}+\Phi_{g})}{\partial
R}\right)+\frac{\partial^{2}(\Phi_{s}+\Phi_{g})}{\partial z^{2}}=4\pi
G(\Sigma_{s}+\Sigma_{g})\delta(z).$ (12)

We now introduce small perturbations in the above basic equations:

$\Sigma_{i}=\Sigma_{0,i}+\epsilon\Sigma_{1,i},\\\
\Phi_{i}=\Phi_{0,i}+\epsilon\Phi_{1,i},\\\ v_{i}=v_{0,i}+\epsilon v_{1,i},\\\
u_{i}=\epsilon u_{1,i}.$ (13)

The quantities $\Sigma_{0,i}$ and $\Phi_{0,i},\,v_{0,i},\,u_{0,i}$ are the
locally unperturbed states, and the perturbed quantities are denoted by
$\Phi_{1,i},\,\Sigma_{1,i},\,v_{1,i},\,u_{1,i}$, where the value of
$\epsilon<<1$. Substituting equation (13) in equations [(9), (10), (11), (12)]
and keeping only the first order terms $[\epsilon^{1}]$, we obtain the
governing equations for the perturbed quantities. But, before that, in order
to better understand how the external potential interacts with the two-
component ${}^{\prime}star+gas^{\prime}$ disc, we write down the zeroth order
terms $[\epsilon^{(0)}]$ corresponding to equation (9):

$\frac{v^{2}_{0,i}}{R}=c^{2}_{i}\frac{\partial\ln\Sigma_{0,i}}{\partial
R}+\frac{\partial\Phi_{0,s}}{\partial R}+\frac{\partial\Phi_{0,g}}{\partial
R}+\frac{\partial\Phi_{ext}}{\partial R}.$ (14)

In the above equation, the contribution of the term
$c^{2}_{i}\partial\ln\Sigma_{0,i}/\partial R$ is negligible since the velocity
dispersion is very small compared to the rotation velocity (Binney & Tremaine,
2011). Further, we write, $v^{2}_{0,s}=R\partial\Phi_{0,s}/\partial R$,
$v^{2}_{0,g}=R\partial\Phi_{0,g}/\partial R$ and
$v^{2}_{ext}=R\partial\Phi_{ext}/\partial R$, we obtain:

$v^{2}_{net}=v^{2}_{0,s}+v^{2}_{0,g}+v^{2}_{ext}.$ (15)

In the above equation, we have labeled $v_{0,i}$ as $v_{net}$ since it
contains the effective contribution from the stars, gas, and the external
potential. The value of $v_{net}$ is typically determined through observations
of neutral hydrogen in galaxies (De Blok et al., 2001; Oh et al., 2011; Lelli
et al., 2016). We express the circular velocity of stars, gas, and the
external potential in terms of the circular frequency as
$v_{0,s}=R\Omega_{0,s}$, $v_{0,g}=R\Omega_{0,g}$ and $v_{ext}=R\Omega_{ext}$.
This leads to
$v^{2}_{net}=R^{2}\Omega^{2}_{net}=R^{2}(\Omega^{2}_{disc}+\Omega_{ext})$,
where $\Omega^{2}_{disc}=\Omega^{2}_{0,s}+\Omega^{2}_{0,g}$.
The stability criterion derived by Toomre (1964), Goldreich & Lynden-Bell
(1965), Jog & Solomon (1984a), Jog (1996) applies exclusively to star and/or
gas disc. The net rotation contains the effect of only star and/or gas but
does not contain the contribution of the external potential to the net
rotation, as shown in equation (15). In their treatment, the centrifugal force
balances the unperturbed potential of either star or/and gas, i.e.,
$v^{2}_{net}/R=\partial(\Phi_{0,s}+\Phi_{0,g})/\partial R$ or
$v^{2}_{net}/R=\partial\Phi_{0}/\partial R$ for a single component, but does
not consider the contribution of the dark matter to the net rotation which
enters our equations as an external potential $(\partial\Phi_{ext}/\partial
R)$. However, when reconstructing stability using observed properties, the
observed rotation curve is used, which includes contributions from stars, gas,
and the dark matter halo. In contrast, the analytic treatment considers
contributions only from the stars and/or gas disk. Jog (2014) identifies this
difference between the analytical treatment and the observational
reconstruction of the stability criterion in the literature and derives a
modified stability criterion for a one-component disc that includes the
contribution of the dark matter halo to the net rotation.
Following the short detour aimed at understanding how the external potential
interacts with the ${}^{\prime}star+gas^{\prime}$ disc, we now write down the
linearized equations governing the growth of perturbed quantities. The first
order terms in $\epsilon^{(1)}$ are given by:

$\frac{\partial u_{1,i}}{\partial
t}-2\Omega_{net}v_{1,i}+\frac{c_{i}^{2}}{\Sigma_{0,i}}\frac{\partial\Sigma_{1,i}}{\partial
R}+\frac{\partial(\Phi_{1,s}+\Phi_{1,g})}{\partial R}=0,$ (16)

$\frac{\partial v_{1,i}}{\partial t}-2B_{net}u_{1,i}=0,$ (17)

$\frac{\partial\Sigma_{1,i}}{\partial t}+\Sigma_{0,i}\frac{\partial
u_{1,i}}{\partial R}=0,$ (18)

and the Poisson equation for the thin disc assumes the form (Toomre, 1964);

$\frac{\partial(\Phi_{1,s}+\Phi_{1,g})}{\partial R}=-2\pi
iG(\Sigma_{1,s}+\Sigma_{1,g}).$ (19)

In equation (16), $v_{net}$ is expressed as $v_{net}=\Omega_{net}R$, and the
term $v^{2}_{net}/R-\partial(\Phi_{0,s}+\Phi_{0,g}+\Phi_{ext})/\partial R$
cancels, as it is just the centrifugal term balancing the total unperturbed
potential of the two-component disc and the external potential. In equation
(17), we have expressed $[\Omega_{net}+\partial(\Omega_{net}R)/\partial
R]=-2B_{net}$, where $B_{net}$ is the Oort constant. Further, substituting
$\Omega_{net}=\sqrt{\Omega^{2}_{disc}+\Omega^{2}_{ext}}$ in
$[\Omega_{net}+\partial(\Omega_{net}R)/\partial R]=-2B_{net}$, it is
straightforward111$\kappa^{2}_{net}=-4B_{net}\Omega_{net}=\left(R\frac{d\Omega^{2}_{disc}}{dR}+4\Omega^{2}_{disc}\right)+\left(R\frac{d\Omega^{2}_{ext}}{dR}+4\Omega^{2}_{ext}\right),\\\
\kappa^{2}_{disc}=\left(R\frac{d\Omega^{2}_{disc}}{dR}+4\Omega^{2}_{disc}\right),\,\kappa^{2}_{ext}=\left(R\frac{d\Omega^{2}_{ext}}{dR}+4\Omega^{2}_{ext}\right)$
to show that $\kappa^{2}_{net}=\kappa^{2}_{disc}+\kappa^{2}_{ext}$, where
$\kappa_{net}$ is the net epicyclic frequency defined as
$\kappa^{2}_{net}=-4B_{net}\Omega_{net}$. In equation (18), the term
$1/R[\partial(Ru_{1,i}\Sigma_{0,i})/\partial R]$ is approximated as
$\Sigma_{0,i}\partial u_{1,i}/\partial R$, as $R\Sigma_{0,i}$ will vary
gradually with $R$ when compared with the rapid oscillatory behaviour of
$u_{1,i}$.

##  3 Dispersion relation in the presence of external field

In this section, we will derive the dispersion relation for the two-component
disc in the presence of an external field of the dark matter halo. We will
show that the linearized equations [(16), (17), (18), (19)] governing the
evolution of the perturbed quantities can be recast to resemble coupled wave
equations with extra terms and thus admit solutions of the form
$e^{ik.r-\omega t}$. Indexing equations [(16), (17), (18), (19)] for stars;

$\frac{\partial u_{1,s}}{\partial
t}-2\Omega_{net}v_{1,s}+\frac{c_{s}^{2}}{\Sigma_{0,s}}\frac{\partial\Sigma_{1,s}}{\partial
R}+\frac{\partial(\Phi_{1,s}+\Phi_{1,g})}{\partial R}=0,$ (20)

$\frac{\partial v_{1,s}}{\partial t}-2B_{net}u_{1,s}=0,$ (21)

$\frac{\partial\Sigma_{1,s}}{\partial t}+\Sigma_{0,s}\frac{\partial
u_{1,s}}{\partial R}=0.$ (22)

Operating with $\partial/\partial R$ on equation (20), and eliminating the
terms $\partial/\partial R\left(\partial u_{1,s}/\partial{t}\right)$ by taking
the time derivative of equation (22), which will give $\partial/\partial
R\left(\partial
u_{1,s}/\partial{t}\right)=\left(-1/\Sigma_{0,s}\right)\partial^{2}\Sigma_{1,s}/\partial
t^{2}$. Similarly, $\partial v_{1,s}/\partial R$ is eliminating by operating
$\partial/\partial R$ on equation (21) and substituting for $\partial
u_{1,s}\partial R$ from equation (22) to get $\partial v_{1,s}/\partial
R=-2B_{net}\Sigma_{1,s}/\Sigma_{0,s}$. And finally substituting for
$\partial^{2}(\Phi_{1,s}+\Phi_{1,g})/\partial R^{2}$ with equation (19), we
obtain:

$\frac{\partial^{2}\Sigma_{1,s}}{\partial
t^{2}}-c_{s}^{2}\frac{\partial^{2}\Sigma_{1,s}}{\partial
R^{2}}-4\Omega_{net}B_{net}\Sigma_{1,s}+2\pi
iG\Sigma_{0,s}\frac{\partial}{\partial R}(\Sigma_{1,s}+\Sigma_{1,g})=0.$ (23)

Similarly, the equation for gas reads ;

$\frac{\partial^{2}\Sigma_{1,g}}{\partial
t^{2}}-c_{g}^{2}\frac{\partial^{2}\Sigma_{1,g}}{\partial
R^{2}}-4\Omega_{net}B_{net}\Sigma_{1,g}+2\pi
iG\Sigma_{0,g}\frac{\partial}{\partial R}(\Sigma_{1,s}+\Sigma_{1,g})=0.$ (24)

The above equations resemble wave equations and will indeed admit plane wave
ansatz. Substituting $e^{i(k.r-\omega t)}$ for the perturbed quantities in
equations (23) and (24), we obtain

$\Sigma_{1,s}=\frac{-2\pi
Gk\Sigma_{0,s}\Sigma_{1,g}}{(\omega^{2}-c_{s}^{2}k^{2}-\kappa^{2}_{net}+2\pi
G\Sigma_{0,s}k)},$ (25)

and similarly

$\Sigma_{1,g}=\frac{-2\pi
Gk\Sigma_{0,g}\Sigma_{1,s}}{(\omega^{2}-c_{g}^{2}k^{2}-\kappa^{2}_{net}+2\pi
G\Sigma_{0,g}k)}.$ (26)

Combining equations (25) and (26), the final dispersion relation reads;

$\displaystyle(\omega^{2}-c_{s}^{2}k^{2}-\kappa^{2}_{net}+2\pi
G\Sigma_{0,s}k)(\omega^{2}-c_{g}^{2}k^{2}-\kappa^{2}_{net}+2\pi
G\Sigma_{0,g}k)=$ (27) $\displaystyle(2\pi G\Sigma_{0,s}k)(2\pi
G\Sigma_{0,g}k).$

By setting the contribution of the external field to zero, i.e.,
$\kappa_{ext}=0$, equation (27) becomes equivalent to the dispersion relation
for a two-component galactic disk, as shown in equation (17) of Jog & Solomon
(1984a). Further, if either of $\Sigma_{0,s},c_{s}=0$ or
$\Sigma_{0,g},c_{g}=0$, equation (27) reduces to the case of a one-component
disc under the influence of the external field (Jog, 2014). In deriving the
above dispersion relation, we have started with a two-component disc in
equilibrium with an external force field of dark matter halo. We then
introduced small perturbations and compose the linearized perturbation
equations, which resemble plane wave equations and then use the plane wave
ansatz to derive the dispersion relation.

##  4 Condition for stability

In this section, we derive the stability criterion for assessing if the two-
component disc in the force field of the dark matter halo is susceptible to
the growth of axisymmetric instabilities or not.
Firstly, we define the following quantities:

$\displaystyle\alpha_{s}=\kappa^{2}_{net}+c_{s}^{2}k^{2}-2\pi G\Sigma_{0,s}k,$
(28) $\displaystyle\alpha_{g}=\kappa^{2}_{net}+c_{g}^{2}k^{2}-2\pi
G\Sigma_{0,g}k,$ $\displaystyle\beta_{s}=2\pi G\Sigma_{0,s}k,$
$\displaystyle\beta_{g}=2\pi G\Sigma_{0,g}k.$

Substituting equation (28) in (27), the dispersion relation and the respective
roots are given by;

$\displaystyle\omega^{4}-\omega^{2}(\alpha_{s}+\alpha_{g})+(\alpha_{s}\alpha_{g}-\beta_{s}\beta_{g})=0$
(29)
$\displaystyle\omega^{2}_{\pm}=\frac{1}{2}(\alpha_{s}+\alpha_{g})\pm\frac{1}{2}((\alpha_{s}+\alpha_{g})^{2}-4(\alpha_{s}\alpha_{g}-\beta_{s}\beta_{g}))^{\frac{1}{2}}.$

For a one-component disc, $\alpha_{g}\geq 0$ or $\alpha_{s}\geq 0$ is the
sufficient condition for stability. For a marginally stable one-component disc
a function $F$ can be defined as $F=2\pi
G\Sigma_{0}k/(\kappa^{2}_{net}+k^{2}c^{2})$. A value of $F=1$ indicates
marginal stability, $F>1$ represents an unstable disc and $F<1$ represents a
stable disc. The value of $k_{min}$ for the one-component disc is obtained by
putting $d\omega^{2}/dk=0$, where $\omega^{2}=\kappa^{2}_{net}-2\pi
G\Sigma_{0}k+c^{2}k^{2}$, which yields $k_{min}=\pi G\Sigma_{0}/c^{2}$.
Evaluating $F$ at $k_{min}$ yields $F=2/(1+Q^{2})$. For the one-component
system in the force field of an external potential, Q is defined as
$Q=\kappa_{net}c/\pi
G\Sigma_{0}=q\sqrt{1+(\kappa^{2}_{ext}/\kappa^{2}_{disc})}$, where
$q=\kappa_{disc}c/\pi G\Sigma_{0}$.
The condition for marginal stability of two-component disc reads
$\omega^{2}_{-}=0$ or $(\alpha_{s}\alpha_{g}-\beta_{s}\beta_{g})=0$, and for
the disc to be unstable the conditions is
$\alpha_{s}\alpha_{g}-\beta_{s}\beta_{g}<0$. With simple algebra the condition
for neutral equilibrium, $(\alpha_{s}\alpha_{g}-\beta_{s}\beta_{g})=0$ can be
written as:

$F=\frac{2\pi
G\Sigma_{0,s}k}{\kappa^{2}_{disc}+\kappa^{2}_{ext}+k^{2}c_{s}^{2}}+\frac{2\pi
G\Sigma_{0,g}k}{\kappa^{2}_{disc}+\kappa^{2}_{ext}+k^{2}c_{g}^{2}}$ (30)

where F=1. In the above, we have expressed
$\kappa^{2}_{net}=\kappa^{2}_{disc}+\kappa^{2}_{ext}$ to gauge the effect of
the external potential on the ${}^{\prime}star+gas^{\prime}$ disc. See, the
discussion following equation (19) for deriving $\kappa_{net}$ in terms of
$\kappa_{disc}$ and $\kappa_{ext}$. We define the gas fraction
$f=\Sigma_{0,g}/(\Sigma_{0,s}+\Sigma_{0,g})$, and
$X_{s-g}=\kappa^{2}_{disc}/[2\pi G(\Sigma_{0,s}+\Sigma_{0,g})k_{min}]$.
$X_{s-g}$ is the dimensionless wavelength at which it is hardest to stabilize
the two-component system. The value of $k_{min}$ for the two-component system
is given by conditions, $d\omega^{2}_{-}/dk=0$, or
$d(\omega^{2}_{+}\omega^{2}_{-})/dk=0$, i.e. finding
$d(\alpha_{s}\alpha_{g}-\beta_{s}\beta_{g})/dk$ which yields;

$\displaystyle k^{3}(4c_{s}^{2}c_{g}^{2})-3k^{2}(2\pi
G\Sigma_{0,s}c_{g}^{2}+2\pi G\Sigma_{0,g}c_{s}^{2})$ (31)
$\displaystyle+2k\kappa^{2}_{net}(c_{g}^{2}+c_{s}^{2})-(2\pi
G\Sigma_{0,s}+2\pi G\Sigma_{0,g})\kappa^{2}_{net}=0$

The function $F$ for the two-component model is a superposition of the one-
component cases (Jog & Solomon, 1984a). Thus, in analogy with the one-
component case, the condition for stability of the two-component disc under
the force field of external potential is defined as:

$\frac{2}{1+Q_{T}^{2}}=\frac{(1-f)}{X_{s-g}(1+\frac{(1-f)^{2}q_{s}^{2}}{4X^{2}_{s-g}}+R)}+\frac{f}{X_{s-g}(1+\frac{f^{2}q_{g}^{2}}{4X^{2}_{s-g}}+R)}.$
(32)

In the above equation, $R$ quantifies the contribution of the external
potential on the two-component ${}^{\prime}star+gas^{\prime}$ and is defined
as $R=\kappa^{2}_{ext}/\kappa^{2}_{disc}$. Also, $q_{s}$ and $q_{g}$ are the
classical one-component stability criterion for stars and gas, defined as
$q_{s}=\kappa_{disc}c_{s}/\pi G\Sigma_{0,s}$ and $q_{g}=\kappa_{disc}c_{g}/\pi
G\Sigma_{0,g}$ respectively. The above condition is equivalent to the
stability condition $Q_{s-g}$ derived by Jog (1996) in the absence of the
external force field $(R=0)$. For the sake of continuity of notation, we
denote the stability criterion for the two-component disc in the absence of an
external field using $q_{T}$. The disc is stable against the growth of
axisymmetric instabilities when $Q_{T}>1$, and the disc is susceptible to the
growth of axisymmetric perturbations when $Q_{T}<1$.

##  5 Results

###  5.1 Marginal stability of one-component disc under the influence of dark
matter halo

To gain better insight into the role of dark matter on a two-component disc,
we first investigate the impact of the force field of dark matter halo in
driving the stability levels in a one-component disc. The dispersion relation
for a one-component disc is given by:

$\omega^{2}=(\kappa^{2}_{disc}+\kappa^{2}_{ext})k^{(0)}+c^{2}k^{(2)}-2\pi
G\Sigma_{0}k^{(1)}.$ (33)

In the above equation, at a large value of $k$, $k^{2}$ will dominate; thus,
pressure stabilizes the disc at small scales. At small $k$, i.e., $k^{(0)}$
the differential rotation of the disc ($\kappa^{2}_{disc}$) and the dark
matter halo ($\kappa^{2}_{ext}$) stabilize the disc at large scales. At
intermediate $k$, the self-gravity of the galactic disc becomes important. The
field due to the external potential $(\kappa^{2}_{ext})$ adds up with the
differential rotation of the disc $(\kappa^{2}_{disc})$ and will stabilize the
disc. Next, we inspect the marginal stability of the one-component galactic
disc. Putting $\omega^{2}=0$, equation (33) can be recast to obtain a
quadratic equation in $k$,

$1+\frac{Q^{2}}{4}\frac{k^{2}}{k^{\prime 2}_{T}}-\frac{k}{k^{\prime}_{T}}=0,$
(34)

where, $Q=q(1+\kappa^{2}_{ext}/\kappa^{2}_{disc})^{\frac{1}{2}}$,
$k^{\prime}_{T}=k_{T}(1+\kappa^{2}_{ext}/\kappa^{2}_{disc})$,
$k_{T}=\kappa^{2}_{disc}/2\pi G\Sigma_{0}$ and defining
$\zeta^{\prime}=k^{\prime}_{T}/k=k_{T}/k(1+\kappa^{2}_{ext}/\kappa^{2}_{disc})$,
i.e $\zeta^{\prime}=\zeta(1+\kappa^{2}_{ext}/\kappa^{2}_{disc})$. With the
above substitutions, equation (34) can be written as
$Q=2\left[\zeta(1+R)\left(1-\zeta(1+R)\right)\right]^{\frac{1}{2}}$, where
$R=\kappa^{2}_{ext}/\kappa^{2}_{disc}$.

Figure 1: The marginal stability of a one-component disc under the influence
of the external force field of dark matter halo.

In Figure. 1, we show the effect of the external force field of the dark
matter halo on the stability of the one-component disc. We find that upon
increasing the contribution of the dark matter by increasing the value of $R$,
the maximum value of $Q$ is shifted towards a smaller value of $\zeta$,
indicating that a larger contribution from dark matter to the total potential
can effectively stabilize the galaxy over large scales. Further, from
$Q=q\sqrt{1+R}$, we can see that when
$R=\kappa^{2}_{ext}/\kappa^{2}_{disc}=0$, the value of $Q$ corresponds to the
classical stability criterion $(q)$ derived by Toomre (1964). The stability
criterion derived by Toomre (1964) considers the self-gravity of only one
component and does not include the contribution of the external potential due
to dark matter halo. The centrifugal force is balanced only by the
corresponding force due to the unperturbed potential of stars/gas. The
marginal stability in the absence of the external potential is given by $q=1$.
The maximum value value of $Q$, when $R=0.5$, or
$\kappa_{ext}=(\kappa_{disc}/\sqrt{2})$ is equal to 1.2 compared to 1 when
R=0, indicating that $Q>q$. Thus, it is evident that the addition of dark
matter to the total potential increases the marginal stability levels and
makes it much harder to destabilize the one-component disc, making the disc
more stable against the growth of instabilities. The stability criterion in
the presence of the external dark matter halo
$Q=q(1+\kappa^{2}_{ext}/\kappa^{2}_{disc})^{\frac{1}{2}}$, can be written as
$Q=\kappa_{net}c/\pi G\Sigma_{0}$. The one-component stability criterion
derived by Toomre (1964) is applicable only for stars/gas. However, we note
that the mathematical expression for the stability criterion in the presence
of an external halo remains unchanged $(Q=\kappa_{net}c/\pi G\Sigma_{0})$.
Thus, when reconstructing $q$ from observations following the classical
treatment by Toomre (1964), the contribution of the external force field of
dark matter is implicitly accounted for. In other words, using net epicyclic
frequency $(\kappa_{net})$ derived from the observed rotation curve in $q$
(Toomre, 1964) is equivalent to computing $Q$ derived in this work.

###  5.2 Role of dark matter on the stability of two-component disc

Figure 2: The effect of varying the external force of dark matter halo on the
stability of the two-component disc. In the top panel, the gas fraction is
fixed at f=0.01, and in the bottom panel, the gas fraction is fixed at f=0.3.
The color bar indicates net stability levels given by $Q_{T}$.

In §5.1, we found that the external potential due to the dark matter halo
increases the marginal stability of the one-component disc, indicating that
effectively $Q>q$, or that it is now harder to destabilize the disc due to the
force field of the dark matter halo. Jog & Solomon (1984a) and Jog (1996),
show that the addition of gas disc makes the galaxy more prone to the growth
of instabilities or in other words, the two-component disc is more unstable
than either component by itself. The stability criterion presented in this
work allows us to assess if gas is enough to lead to the growth of local
instabilities, even in the presence of a stabilizing dark matter component. In
Figure 2, we present the two-component stability criterion $Q_{T}$ as a
function of $q_{s}$ and $q_{g}$, respectively. In the top panel, we have fixed
the value of the gas fraction at $f=0.01$; in the bottom panel, we have fixed
the gas fraction at 0.3. When the external force field due to the dark matter
halo is zero $R=0$, we find that at a fixed value of the gas fraction, the
value of the $q_{T}$ is lower than the values of $q_{s}$ and $q_{g}$. This
supports the earlier findings by Jog & Solomon (1984a); Jog (1996), which show
that a two-component disc is less stable than a disc composed only of stars or
gas. The results indicate that the two-component disc is more prone to the
development of gravitational instabilities than a single-component disc. For
example, when $f=0.01,\,R=0$, value of $q_{T}=1.95$, when $q_{s}$ and $q_{g}$
are equal to 2.5. Similarly, when we increase the gas fraction to $f=0.3$,
keeping $R=0$, $q_{T}$ now becomes 1.35 when $q_{s},q_{g}=2.5$. Further, when
$f=0.3\,,R=0$, and $q_{s}\,,q_{g}=1.5$, the value of $q_{T}$ drops to 0.75.
This shows that adding a second disc in the absence of external potential due
to dark matter effectively renders the two-component disc susceptible to the
growth of gravitational instabilities, even though the stars and gas are
stable by themselves.
We will now discuss the effect of dark matter on the stability of the two-
component disc by varying the value of $R$. It is evident from Figure 2 that
for a given gas fraction, when we move from left to right, the value of the
marginal stability of the two-component disc increases with increasing R. For
example, at a gas fraction equal to 0.01, when both $q_{s}$ and $q_{g}$ are
equal to 2.5, $q_{T}$ is equal to 1.95 for $R=0$. However, upon increasing the
contribution of dark matter to the total potential (i.e., $R=1$), $Q_{T}$
becomes 2.4. A similar effect is observed at a higher gas fraction, when
$f=0.3$ and both $q_{s}$ and $q_{g}$ are equal to 2.5, the value of $q_{T}$ is
1.35 in the absence of dark matter ($R=0$). However, when the contribution of
dark matter is included ($R=1$), $Q_{T}=1.8$. The external potential of the
dark matter halo stabilizes the two-component system, which would otherwise be
prone to the growth of axisymmetric instabilities. At smaller values of
$q_{s}$ and $q_{g}$ equal to 1.5 and a gas fraction equal to 0.3, the two-
component disc becomes susceptible to the growth of gravitational
instabilities ($q_{T}=0.75$) when $R=0$. However, upon including the
contribution of the dark matter halo $(R=1)$ at $f=0.3$, the two-component
system stabilizes itself $(Q_{T}>1)$. Thus, we note that an external force due
to the dark matter halo effectively suppresses the growth of local
axisymmetric instabilities. However, the two-component system can be
susceptible to axisymmetric instabilities on rarer occasions, even in the
presence of stabilizing external potential due to the dark matter halo. An
example is provided in Figure 2, when $R=0.5$, and the gas fraction is equal
to $0.3$ and $q_{s},\,q_{g}\leq 1.5$. The two-component system has $Q_{T}\leq
1$, indicating that the system is prone to growth of local axisymmetric
instabilities.

##  6 Application

Figure 3: In the first row, we show stability analysis for model of Milky Way.
In the second row, we show the stability analysis for a model of a low-mass
stellar disc embedded in a massive dark matter halo. In the third row, we
depict the analysis for a low-mass stellar disc embedded in a low-mass dark
matter halo akin to a low surface brightness galaxy. The fourth panel shows
the stability analysis for a baryon-dominated galaxy akin to cold rotating
disc galaxies observed in the early universe. In realistic scenarios, $Q_{g}$,
$Q_{s}$, and $Q_{T}$ correspond to the stability curves estimated from
observations. The red dashed line indicates marginal stability levels.

From the analysis presented in the previous section, we understand that the
external force field of the dark matter halo stabilizes the two-component
system of $stars+gas$. However, in rare instances where the force due to the
dark matter halo is insufficient compared to the destabilizing effect of the
gas disc, the two-component system may become prone to the growth of local
axisymmetric instabilities. In this section, we will investigate the role of
dark matter on the stability of two-component models of nearby galaxies like
the Milky Way and low surface brightness galaxies, and models of galaxies
observed in the early universe. The stellar distribution in our galaxy follows
an exponential surface density given by

$\Sigma_{s}(R)=\Sigma_{s0}e^{-R/R_{D}},$ (35)

where $\Sigma_{s0}=640M_{\odot}pc^{-2}$ and $R_{D}=3.2kpc$ (Mera et al., 1998)
are the central surface density and the disc scalelength. The gas distribution
in the Galaxy is given by

$\Sigma_{g}(R)=\Sigma_{g0}e^{-1.65R/R_{25}},$ (36)

in the above equation $\Sigma_{g0}=28.2M_{\odot}pc^{-2}$ is the central
density of the gas disc and $R_{25}$ is the radius at which the B-band surface
brightness drops to $25.5\,mag\,arsec^{-2}$, $R_{25}=4R_{D}$ (Bigiel & Blitz,
2012). The stellar velocity dispersion is given by (Leroy et al., 2008; Romeo
& Mogotsi, 2017) $\sigma_{s}(R)=(1/0.6)\sqrt{(2\pi GR_{D}\Sigma_{s}(R)/7.3}$,
and we use a constant gas velocity dispersion equal to $10kms^{-1}$ (Tamburro
et al., 2009; Mogotsi et al., 2016). The circular velocity $(v_{c})$
corresponding to the exponential distribution is given as (Binney & Tremaine,
2011)

$v^{2}_{c}(R)=4\pi G\Sigma_{0}R_{D}y^{2}[I_{0}(y)K_{0}(y)-I_{1}(y)K_{1}(y)],$
(37)

where $y=R/2R_{D}$ and $I_{0},I_{1}$ and $K_{0},K_{1}$ are the modified Bessel
functions of the first and second kind. The epicyclic frequency $\kappa$ at a
radius R is defined as

$\kappa^{2}(R)=\bigg{(}R\frac{d\Omega^{2}(R)}{dR}+4\Omega^{2}(R)\bigg{)},$
(38)

where $\Omega$ is the angular frequency defined as
$\Omega^{2}(R)=\frac{v^{2}_{c}}{R^{2}}$. The dark matter density is given by a
pseudo-isothermal halo, with a central density
$\rho_{0}=0.035M_{\odot}pc^{-3}$ and a core radius $R_{c}=5kpc$ (Mera et al.,
1998). The epicyclic frequency due to the pseudo-isothermal dark matter halo
is given by

$\kappa^{2}_{PIS}(R)=4\pi
G\rho_{0}\left[\frac{2R_{c}^{2}}{R^{2}+R_{c}^{2}}+\frac{R_{c}^{4}}{R^{2}(R^{2}+R_{c}^{2})}-\frac{R_{c}^{3}}{R^{3}}\tan^{-1}\left(\frac{R}{R_{c}}\right)\right].$
(39)

With all the building blocks needed to compute the stability in place, we will
now discuss the role of dark matter halo on different galaxy models.
Case 1: Stability of Milky Way
We show the stability analysis of the Milky Way in the first row of Figure 3.
We find the minimum value of the stability for stars and gas is
$q^{min}_{s}=1.5$ and $q^{min}_{g}=1.8$, indicating that individually stars
and gas are stable. However, the two-component formalism yields a
$q^{min}_{T}=0.88$, indicating that the two-component system is prone to the
growth of local gravitational instabilities in the absence of the external
potential due to the dark matter halo. Now, the addition of the dark matter to
the total potential naturally increases the one-component stability from
$q^{min}_{s}=1.5$ to $Q^{min}_{s}=2.7$ and $q^{min}_{g}=1.8$ to
$Q^{min}_{g}=3.8$. Moreover, the two-component system, which was unstable with
$q^{min}_{T}=0.88$, now has $Q^{min}_{T}=1.7$, highlighting the importance of
the external potential of the dark matter halo in stabilizing massive disc
galaxy like Milky Way.
Case 2: Low mass disc in a high mass halo
In order to better gauge the role of dark matter on the stability levels of
the two-component disc, we minimize the contribution of the stars and gas disc
to the total potential. We lower the stellar and the gas surface density to
typical values observed in the low surface brightness galaxies;
$\Sigma_{s0}=100M_{\odot}pc^{-2}$ and $R_{D}=2.5kpc$ and the gas surface
density to $14.2M_{\odot}pc^{-2}$ (De Blok et al., 2001; Di Paolo et al.,
2019). We keep the values of the dark matter halo to that of the Milky Way. We
show the rotation velocity for this mass distribution in the second row of
Figure 3. We can see that $\kappa_{ext}$ is significantly higher than
$\kappa_{disc}$, highlighting that the dark matter is the dominant mass
component. We find that $q^{min}_{s}=1.7$ and $q^{min}_{g}=1.8$, indicating
that the stars and gas are stable on their own. However, similar to the Milky
Way, the two-component star+gas system has $q^{min}_{T}=0.9$, making the disc
susceptible to the growth of local gravitational instabilities in the absence
of the dark matter halo. However, upon including the contribution of the dark
matter halo, we find $Q^{min}_{s}=6.9$, $Q^{min}_{g}=8$, and
$Q^{min}_{T}=4.4$, indicating that a higher contribution of the dark matter to
the total potential is reflected in higher net stability levels of the two-
component system.
Case 3: Low mass disc in a low mass halo
As a final example, we will inspect the effect of a low-mass stellar and gas
disc embedded in a low-mass dark matter halo akin to the mass distribution of
a low surface brightness galaxy. We keep the surface density of the stars,
gas, and dark matter halo parameters to the typical values obtained from the
mass models of low surface brightness galaxies:
$\rho_{0}=0.066M_{\odot}pc^{-3}$ and $R_{c}=1.5kpc$ (De Naray et al., 2008; Di
Paolo et al., 2019). The parameters for the stars and gas are the same as in
Case 2. We show the rotation velocity corresponding to this mass distribution
in the third row of Figure 3. The minimum value of $q_{s}$, $q_{g}$ and
$q_{T}$ are comparable to values obtained in case 2: $q^{min}_{s}=1.7$,
$q_{g}=1.8$ and $q_{T}=0.9$. However, since the contribution of dark matter to
the total potential is small compared to the massive dark matter halo of the
Milky Way, the shift in the stability curves upon adding a dark matter halo is
also small. We find that $Q^{min}_{s}=3.2$, $Q^{min}_{g}=3.6$ and
$Q^{min}_{T}=2.1$. In both Case 2 and Case 3, we find that the force due to
the dark matter potential stabilizes the two-component low surface brightness
disc, which is otherwise unstable. The only difference is that a massive halo
contributes significantly to the overall stability. This aligns with the
previous finding in Garg & Banerjee (2017); Aditya (2023), which shows that
dark matter is important in regulating the stability of low surface brightness
galaxies.

##  7 Discussion

A large number of recent studies show that the galaxies observed at high
redshift are dominated by baryons (Rizzo et al., 2021; Genzel et al., 2020,
2017, 2014). We construct a galaxy model in which the contribution of stars
and gas exceeds that of the dark matter halo in the total mass budget. The
rotation curve decomposition by Rizzo et al. (2021) show that stellar disc
makes maximum contribution to the total rotation curve. However, the
contribution from the dark matter halo and the gas disc are typically
comparable. In our model, the stellar disc has a surface density profile
comparable to the Milky Way: $\Sigma_{0}=640M_{\odot}pc^{-2}$ and
$R_{d}=3.2kpc$. We increase the gas surface density from
$28.2M_{\odot}pc^{-2}$ for the Milky Way to $200M_{\odot}pc^{-2}$ and keep the
scalelength comparable to the stellar disc. However, the gas disc continues to
be a cold component with a velocity dispersion of $10km/s$. We also lower the
contribution of the dark matter by reducing the dark matter density and core
radius to $0.05M_{\odot}pc^{-3}$ and $2kpc$, respectively. We aim to ascertain
the contribution of the dark matter halo to the net stability levels in
baryon-dominated systems akin to the cold rotating disc galaxies observed in
the early universe. We show the results in the fourth panel of Figure 3. We
find that in the absence of potential due to the dark matter halo,
$q^{min}_{s}=1.5$, $q^{min}_{g}=0.8$ and $q^{min}_{T}=0.5$, indicating that a
massive cold two-component disc is susceptible to the growth of local
gravitational instabilities. Although the dark matter halo increases the net
stability levels, the two-component system is still susceptible to local
gravitational instabilities (Bacchini et al., 2024; Aditya, 2023),
$Q^{min}_{s}=2.8$, $Q^{min}_{g}=0.9$, and $Q^{min}_{T}=0.8$. This indicates
that despite the stabilizing nature of the dark matter halo, the net
contribution of the dark matter is insufficient to stabilize the baryon-
dominated cold disc galaxies in the early universe.

##  8 Conclusions

In this study, we have derived detailed theoretical formalism to understand
the role of dark matter and gas fraction on the stability of the two-component
model of galactic disc. We model the galaxy as a coplanar and a coaxial system
of stars and gas in equilibrium with an external dark matter halo. We derive
the equations governing the growth rate of perturbation and, finally, present
a simple stability criterion for appraising the stability of the two-component
disc under the influence of dark matter halo. We find that:

  1. 1.

The two-component disc is more susceptible to the growth of gravitational
instabilities than the individual components. Increasing the gas fraction at a
fixed value of external potential lowers the stability of the two-component
disc, highlighting the role of cold gas in destabilizing the galaxy consistent
with the earlier finding of Jog & Solomon (1984a).

  2. 2.

The external field due to the dark matter halo acts as a stabilizing agent and
increases the net stability levels of the two-component system. In dark
matter-dominated systems, the gravitational force exerted by the dark matter
halo stabilizes the two-component system, even when the system is locally
unstable (Jog, 2014). This indicates that the cold gas component cannot
destabilize the two-component disc when the dark matter halo dominates the
mass budget of the galaxies.

  3. 3.

We apply the stability criterion to the models of the Milky Way and low
surface brightness galaxies and find that the Milky Way and the low surface
brightness discs are locally unstable, when contribution of dark matter is not
included in the total potential. However, the addition of the dark matter to
the total potential significantly increases the net stability levels in these
galaxies (Aditya, 2023). We note that when the contribution of dark matter to
the total mass budget is small, the corresponding effect on the net stability
levels would also be diminished.

  4. 4.

In rare cases, the two-component system can be susceptible to the growth of
gravitational instabilities despite the presence of a stabilizing dark matter
halo potential. One example is found in baryon-dominated cold rotating disc
galaxies observed in the early universe. The influence of dark matter on the
overall gravitational potential is insufficient to stabilize the galaxies
observed in early universe.

##  9 acknowledgements

Aditya would like to thank the referee for their insightful comments that
improved the quality of this manuscript.

##  10 Data Availability

No new data was generated in this work.

## References

  * Aditya (2023) Aditya K., 2023, Monthly Notices of the Royal Astronomical Society, 522, 2543
  * Agertz et al. (2015) Agertz O., Romeo A. B., Grisdale K., 2015, Monthly Notices of the Royal Astronomical Society, 449, 2156
  * Bacchini et al. (2024) Bacchini C., et al., 2024, arXiv preprint arXiv:2405.00103
  * Bigiel & Blitz (2012) Bigiel F., Blitz L., 2012, The Astrophysical Journal, 756, 183
  * Binney & Tremaine (2011) Binney J., Tremaine S., 2011, Galactic dynamics. Vol. 20, Princeton university press
  * De Blok et al. (2001) De Blok W., McGaugh S. S., Rubin V. C., 2001, The Astronomical Journal, 122, 2396
  * De Naray et al. (2008) De Naray R. K., McGaugh S. S., De Blok W., 2008, The Astrophysical Journal, 676, 920
  * Di Paolo et al. (2019) Di Paolo C., Salucci P., Erkurt A., 2019, Monthly Notices of the Royal Astronomical Society, 490, 5451
  * Dopita & Ryder (1994) Dopita M. A., Ryder S. D., 1994, The Astrophysical Journal, 430, 163
  * Elmegreen (1995) Elmegreen B., 1995, Monthly Notices of the Royal Astronomical Society, 275, 944
  * Forbes et al. (2014) Forbes J. C., Krumholz M. R., Burkert A., Dekel A., 2014, Monthly Notices of the Royal Astronomical Society, 438, 1552
  * Garg & Banerjee (2017) Garg P., Banerjee A., 2017, Monthly Notices of the Royal Astronomical Society, 472, 166
  * Genzel et al. (2014) Genzel R., et al., 2014, The Astrophysical Journal, 785, 75
  * Genzel et al. (2017) Genzel R., et al., 2017, Nature, 543, 397
  * Genzel et al. (2020) Genzel R., et al., 2020, The Astrophysical Journal, 902, 98
  * Goldreich & Lynden-Bell (1965) Goldreich P., Lynden-Bell D., 1965, Monthly Notices of the Royal Astronomical Society, 130, 125
  * Hoffmann & Romeo (2012) Hoffmann V., Romeo A. B., 2012, Monthly Notices of the Royal Astronomical Society, 425, 1511
  * Iye (1978) Iye M., 1978, Publications of the Astronomical Society of Japan, 30, 223
  * Jog (1996) Jog C. J., 1996, Monthly Notices of the Royal Astronomical Society, 278, 209
  * Jog (2014) Jog C. J., 2014, The Astronomical Journal, 147, 132
  * Jog & Solomon (1984b) Jog C. J., Solomon P. M., 1984b, The Astrophysical Journal, 276, 127
  * Jog & Solomon (1984a) Jog C., Solomon P., 1984a, The Astrophysical Journal, 276, 127
  * Kalnajs (1983) Kalnajs A. J., 1983, in Symposium-International Astronomical Union. pp 109–116
  * Krumholz & Burkert (2010) Krumholz M., Burkert A., 2010, The Astrophysical Journal, 724, 895
  * Lelli et al. (2016) Lelli F., McGaugh S. S., Schombert J. M., 2016, The Astronomical Journal, 152, 157
  * Leroy et al. (2008) Leroy A. K., Walter F., Brinks E., Bigiel F., De Blok W., Madore B., Thornley M., 2008, The astronomical journal, 136, 2782
  * Lin & Shu (1987) Lin C., Shu F. H., 1987, in , Selected Papers of CC Lin with Commentary: Vol. 1: Fluid Mechanics Vol. 2: Astrophysics. World Scientific, pp 561–570
  * Meidt (2022) Meidt S. E., 2022, The Astrophysical Journal, 937, 88
  * Mera et al. (1998) Mera D., Chabrier G., Schaeffer R., 1998, Astronomy and Astrophysics, 330, 953
  * Mogotsi et al. (2016) Mogotsi K., de Blok W., Caldú-Primo A., Walter F., Ianjamasimanana R., Leroy A., 2016, The Astronomical Journal, 151, 15
  * Oh et al. (2011) Oh S.-H., De Blok W., Brinks E., Walter F., Kennicutt R. C., 2011, The Astronomical Journal, 141, 193
  * Pandey & Van De Bruck (1999) Pandey U., Van De Bruck C., 1999, Monthly Notices of the Royal Astronomical Society, 306, 181
  * Rafikov (2001) Rafikov R. R., 2001, Monthly Notices of the Royal Astronomical Society, 323, 445
  * Rizzo et al. (2021) Rizzo F., Vegetti S., Fraternali F., Stacey H. R., Powell D., 2021, Monthly Notices of the Royal Astronomical Society, 507, 3952
  * Romeo & Falstad (2013) Romeo A. B., Falstad N., 2013, Monthly Notices of the Royal Astronomical Society, 433, 1389
  * Romeo & Fathi (2016) Romeo A. B., Fathi K., 2016, Monthly Notices of the Royal Astronomical Society, 460, 2360
  * Romeo & Mogotsi (2017) Romeo A. B., Mogotsi K. M., 2017, Monthly Notices of the Royal Astronomical Society, 469, 286
  * Romeo & Wiegert (2011) Romeo A. B., Wiegert J., 2011, Monthly Notices of the Royal Astronomical Society, 416, 1191
  * Sellwood (2000) Sellwood J., 2000, Astrophysics and Space Science, 272, 31
  * Sellwood (2012) Sellwood J., 2012, The Astrophysical Journal, 751, 44
  * Tamburro et al. (2009) Tamburro D., Rix H.-W., Leroy A., Mac Low M.-M., Walter F., Kennicutt R., Brinks E., De Blok W., 2009, The Astronomical Journal, 137, 4424
  * Toomre (1964) Toomre A., 1964, The Astrophysical Journal, 139, 1217
  * Toomre (1977) Toomre A., 1977, Annual Review of Astronomy and Astrophysics, 15, 437
  * Wang & Silk (1994) Wang B., Silk J., 1994, The Astrophysical Journal, 427, 759


# Time-Varying Wireless Power Transfer Systems for Improving Efficiency

X. Wang1,2 I. Krois3 N. Ha-Van1 M. S. Mirmoosa1 P. Jayathurathnage1,4 S.
Hrabar3 S. A. Tretyakov1 1Department of Electronics and Nanoengineering, Aalto
University, Finland  2Institute of Nanotechnology, Karlsruhe Institute of
Technology, Germany  3Faculty of Electrical Engineering and Computing,
University of Zagreb, Croatia  4Danfoss Drives, Finland

###### Abstract

Conventional wireless power transfer systems are linear and time-invariant,
which sets fundamental limitations on their performance, including a tradeoff
between transfer efficiency and the level of transferred power. In this paper,
we introduce and study a possibility of temporal modulation for inductive
wireless power transfer systems and uncover that this tradeoff is avoided as a
consequence of varying the inductive coupling strength in time. Our
theoretical analysis reveals that under the optimal modulation depth and
phase, the time modulation can yield a substantial improvement in the WPT
efficiency, while the received power at the load is also improved compared to
the static WPT reference system. We experimentally demonstrate the concept
with a low-frequency system and observe a threefold improvement in efficiency
over the reference static counterpart. This technical capability reconciles
the inherent tradeoff between the WPT efficiency and transferred power, paving
the way for simultaneous advancements in both efficiency and delivered power.

##  I Introduction

Wireless power transfer (WPT) is a promising technology that enables
transmission of electrical power without the need for wires or cables [1, 2,
3, 4]. It is desired that both the transfer efficiency and the transferred
power are high. However, one of the major challenges of WPT designs is a
relatively low efficiency, which can result in significant power losses during
transmission. Moreover, there is always a tradeoff between efficiency and
transferred power [5, 6].

Typically, the WPT efficiency is highly influenced by two critical factors:
the mutual coupling coefficient $k$ and the quality factor $Q$ of the coils
[7, 8, 9]. The product $kQ$ plays a crucial role as it directly characterizes
power losses and, consequently, the overall power transfer efficiency of the
system [10]. Researchers have made efforts to enhance efficiency by developing
various methods to increase these two factors. Recent studies highlighted the
use of multiple-coil WPT systems, such as four-coil systems [2, 11], or relay
resonator wireless systems [5], which leverage the advantages of cascading
multiple mutually coupled coils. For example, in [11], the four-coil system
provides extra mutual coupling to compensate for the efficiency loss when the
transmitting and receiving modules are weakly coupled or placed far apart.
Compared to basic two-coil systems, these multiple-coil systems offer
substantial improvement in system efficiency. However, the matching impedance
condition in these systems still poses an inherent limitation, with the system
energy efficiency never exceeding 50%, which is known as a fundamental limit
of WPT systems. In addition, the coupling coefficient, which is directly
proportional to the mutual inductance, is mainly dependent on the geometry of
the coil structure and the distance between the transmitter and receiver. As a
result, there is a limited room for improvement in this critical factor. In
loosely coupled systems, when the coupling coefficient $k$ is low, it is
advisable to use windings with a high $Q$ factor to increase the power
transfer efficiency of the WPT systems. In other words, winding technology
with low resistance at high-frequency operation is needed for efficient WPT
[12]. Usually, Litz wire consisting of many thin wire strands is used to
reduce the AC winding resistance and improve the $Q$ factor, and therefore
enhance the power transfer efficiency [13, 14].

Nevertheless, conventional WPT systems face inherent limitations in terms of
power transfer and efficiency, particularly when the coupling coefficient and
quality factor are highly restricted by predetermined working conditions.
Typically, beyond a certain threshold of efficiency improvement, there is a
diminishing return in transferred power, contradicting the very purpose of
enhancing efficiency [5, 15]. Balancing these two metrics has become the main
challenge in designing wireless power transfer systems for specific
applications. For example, in some cases, such as in medical implants or
electronic devices [16], maximizing efficiency may be more important than
delivering a high amount of power. In other cases, like electric vehicle
charging, a balance between efficiency and power delivery is sought to ensure
practical usability [17]. Many studies have aimed to enhance the output power
while maintaining the efficiency. This has been achieved through adaptive
adjustments of the operating or resonant frequencies [18, 19]. However,
dynamically adjusting such a system requires additional information about the
operating environment and the efficiency can not be higher than $50\%$.
Essentially, for a fixed WPT system, there is no available method to overcome
these fundamental drawbacks of WPT systems.

Since recently, an alternative way to control the response of electromagnetic
systems by varying the effective properties of those systems in time is
quickly developing [20, 21]. New possibilities arise because modulation of the
system parameters gives an additional degree of freedom that can be used for
enhancing functionalities and breaking conventional limitations [22]. For
example, nonreciprocity [23], frequency translation [24], (parametric)
amplification [25], wideband impedance matching and increasing bandwidth [26,
27] indicate only a small part of this enticing and emerging research
direction whose roots come back to the previous century (see, e.g., Ref.
[21]).

We expect that by incorporating the concept of time modulation as an
additional degree of freedom, it is possible to break the fundamental
limitations of WPT systems. In our previous work [28], we have predicted that
temporally modulating mutual inductance can reconcile the tradeoff between
transferred power and efficiency. In this study, we take a closer look at
interactions of time-varying coupling coils, investigating the roles of not
only the mutual inductance variations, but all the time-varying inductive
components. Based on our theoretical findings, we make a proof-of-concept
experiment to verify the theoretical results.

The paper is organized as follows: In Section II, we present a developed
analytical method that accounts for time-modulated inductance and mutual
inductance at the time-varying circuit level. In Section III, we investigate
the effects of each inductive circuit component of the set of inductively
coupled coils, elucidating the physical meanings of time-varying components
that contribute to possible enhancements of the WPT performance. In Section
IV, we report on a test experimental implementation of a time-modulated WPT
system. A simple WPT configuration based on a low-frequency mechanical
transducer was implemented as a proof-of-concept of time-varying mutual
inductance enhancing WPT efficiency.

##  II Theory

###  II.1 Performance indicators of classical wireless power transfer systems

To establish a basis for comparison, we initially provide a concise overview
of classical WPT systems without any time-modulated components. Figure 1 shows
the circuit diagram of a classical inductively coupled WPT system. The system
comprises a harmonic voltage source $V_{\rm s}(t)$ with internal impedance
$Z_{\rm s}$, a transmitting resonator (depicted as an $LCR$ circuit $L_{1}$,
$C_{1}$, and $R_{1}$), a receiving resonator ($L_{2}$, $C_{2}$, and $R_{2}$),
and a resistive load $R_{L}$. The transmitting resonator is inductively
coupled to the receiving resonator, with the level of coupling quantified by
the mutual inductance ${M}$. First, we assume that all circuit components
remain static, that is, $L_{1,2}(t)$ and $M(t)$ do not depend on time, as in
all conventional WPT systems.

The key performance indicators are the delivered amount of power and the power
transfer efficiency, both of which strongly depend on the values of ${M}$ and
$R_{\rm load}$. At the resonance frequency $\omega_{\rm s}$, the input
impedance $Z_{\rm in}$ seen by the source is given by $Z_{\rm
in}=R_{1}+(\omega_{\rm s}{M})^{2}/(R_{2}+R_{L})$ [29]. The power transferred
from the source to the load $P_{L}$ is maximized when $Z_{\rm in}=Z_{\rm
s}^{*}$ (where $*$ denotes complex conjugate), corresponding to a specific
value of ${M}$ for a given $R_{\rm load}$. Deviations from this mutual
inductance value result in suboptimal power transfer due to impedance mismatch
at the source terminal, leading to reduced delivered power.

Figure 1: Schematic view of an inductively coupled WPT system where all the
inductive components are modulated in time.

The power transfer efficiency $\eta$, primarily influenced by dissipation in
the lossy components and the power source. It can be defined as [5]

$\eta=\frac{P_{L}}{P_{\rm
in}}=\dfrac{P_{L}}{P_{{R1}}+{P_{R2}+{P_{s}}+{P_{L}}}},$ (1)

where $P_{R1}$, $P_{R2}$, and $P_{s}$ represent the power losses in the coil
resistances $(R_{1},R_{2})$ and the source resistance $({\rm{Re}}[Z_{\rm
s}])$. Note that the sum of all the losses and the load power equals to the
total input power to the system.

###  II.2 Time-modulated inductor and mutual inductors in a circuit

The fundamental trade-off between maximizing delivered power and efficiency
[5] follows from passivity and time-invariance of the system, meaning that
power is supplied only by the primary source and all the circuit components do
not change in time. However, if some circuit components are active and/or
change in time under action of some external force, this classical limitation
can be overcome. Here, we focus on the use of time-varying (parametric)
systems. Before we present the results, in this section, we outline the
analytical method to analyze time-varying circuits using the mode-matching
method.

In what follows, we discuss how time-varying circuit components can be
characterized and how currents and voltages in time-varying circuits can be
calculated, with the focus on a time-modulated inductor [28]. Let us consider
an inductor whose inductance $L$ is periodically modulated in time. We assume
that the modulation of inductance is an arbitrary time-periodic function,
$L(t)=L(t+T)$, where $T$ is the modulation period. An arbitrary periodic
function can be expanded into Fourier series as

$L(t)=\sum_{m=-\infty}^{+\infty}l_{m}e^{jm\omega_{\rm M}t},$ (2)

where $\omega_{\rm M}=2\pi/T$ is the fundamental angular frequency of
modulation. Under a time-harmonic excitation, the time-varying inductance is a
parametric device, where an infinite number of frequency harmonics is created.
The frequency of harmonic $n$ is denoted as $\omega_{n}=\omega_{\rm
s}+n\omega_{\rm M}$ where $\omega_{\rm s}$ is the source frequency. Therefore,
the current through a time-varying inductor and voltage across it in time
domain can be generally written as sums of all the frequency components:

$I(t)=\sum_{n=-\infty}^{+\infty}i_{n}e^{j\omega_{n}t},$ (3a)
$V(t)=\sum_{n=-\infty}^{+\infty}v_{n}e^{j\omega_{n}t},$ (3b)

where $i_{n}$ and $v_{n}$ are the complex amplitudes of the current and
voltage harmonics, respectively. We omit the operator Re (real part) in front
of the summation sign, for mathematical convenience in the following
derivations. The voltage across an inductor is the time derivative of its
magnetic flux:

$V(t)=\frac{d\Psi(t)}{dt}=\frac{d[L(t)I(t)]}{dt},$ (4)

where $\Psi$ is the magnetic flux through the inductor. The above differential
equation can be written in integral form after integrating both sides with
respect to time:

$\int V(t)dt=L(t)I(t).$ (5)

After substituting Eqs. (2), (3a), and (3b) into Eq. (5) we have

$\sum_{n=-\infty}^{+\infty}\frac{v_{n}}{j\omega_{n}}e^{j\omega_{n}t}=\sum_{n=-\infty}^{+\infty}\sum_{m=-\infty}^{+\infty}i_{n}l_{m}e^{j\omega_{n+m}t}.$
(6)

In the right side of Eq. (6) we can shift the summation index $n$ by $m$,
i.e., replace $n$ by $n-m$, so that the left and right sides share the same
basis $e^{j\omega_{n}t}$. After removing the summation over $n$ at both sides,
the current-voltage relation for harmonic $n$ is found in the form

$v_{n}=\sum_{m=-\infty}^{+\infty}j\omega_{n}i_{n-m}l_{m}.$ (7)

We can see that due to modulation, the $n$-th order voltage harmonic is a
weighted sum of all the current harmonics $n-m$, where
$m\in[-\infty,+\infty]$, with the impedance coefficients $j\omega_{n}l_{m}$.
If we consider a finite number of harmonics from $n=-N$ to $n=+N$ ($N$ is an
integer), we can write $2N+1$ such linear equations as Eq. (7), and this set
of linear equations can be set in matrix form

$\mathbf{v}=\bar{\bar{Z}}_{L}\cdot\mathbf{i},$ (8)

where

$\bar{\bar{Z}}_{L}=\begin{pmatrix}j\omega_{-N}l_{0}&j\omega_{-N}l_{-1}&\cdots&j\omega_{-N}l_{-2N}\\\
j\omega_{{-N+1}}l_{1}&j\omega_{{-N+1}}l_{0}&\cdots&j\omega_{{-N+1}}l_{1-2N}\\\
\vdots&\vdots&\ddots&\vdots\\\
j\omega_{N}l_{2N}&j\omega_{N}l_{2N-1}&\cdots&j\omega_{N}l_{0}\end{pmatrix}$
(9)

is the impedance matrix of the time-varying inductance which should be
invertible (the inverse of an impedance matrix is an admittance matrix), and

$\displaystyle\mathbf{v}=$
$\displaystyle[v_{-N},v_{-N+1},\cdots,v_{0},\cdots,v_{N-1},v_{N}]^{T}$ (10a)
$\displaystyle\mathbf{i}=$
$\displaystyle[v_{-N},v_{-N+1},\cdots,v_{0},\cdots,v_{N-1},v_{N}]^{T}$ (10b)

denote the voltage and current arrays, respectively. We can see that in time-
modulated circuits the currents and voltages are characterized by arrays
containing complex amplitudes of all the considered frequency components, and
the impedances are matrices. Similarly, a time-varying mutual inductance can
be represented as such matrix, which we denote as $\bar{\bar{Z}}_{M}$

In such multi-frequency systems, the impedances of static components should
also be characterized by matrices. In contrast to time-varying inductances,
static components do not create mutual coupling between different harmonics of
voltages and currents. Therefore, their impedance matrix representations are
purely diagonal. The diagonal terms are the impedance values at the
frequencies of the corresponding harmonic components. We denote the impedance
matrices of components $Z_{s}$, $C_{1}$, $R_{1}$, $C_{2}$, $R_{2}$, and
$R_{L}$ as $\bar{\bar{Z}}_{s}$, $\bar{\bar{Z}}_{C1}$, $\bar{\bar{Z}}_{R1}$,
$\bar{\bar{Z}}_{C2}$, $\bar{\bar{Z}}_{R2}$, and $\bar{\bar{Z}}_{RL}$.

Figure 2: Equivalence of time-varying magnetically coupled coils to time-
varying T-circuits.

Next, we calculate the current and voltage of the circuit considering all the
harmonics. Before doing this, we transform the inductive coupling circuit into
an equivalent T-circuit, as shown in Fig. 2. Such transformation is usually
performed for static systems. In Appendix, we prove that this equivalence is
also valid if the coupling transformer is varying in time ($L_{1,2}(t)$ and
$M(t)$). In this way, the calculation of currents and voltages through or
across each component can follow the traditional way, simply replacing the
scalar parameter values with the corresponding matrices. First, we calculate
the current in the transmitter coil, $\mathbf{i}_{s}$. According to the
equivalent T-circuit in Fig. 2, the total impedance seen from the source can
be calculated as

$\displaystyle\bar{\bar{Z}}_{\rm tot}$
$\displaystyle=\bar{\bar{Z}}_{s}+\bar{\bar{Z}}_{R1}+\bar{\bar{Z}}_{C1}+\bar{\bar{Z}}_{L1}-\bar{\bar{Z}}_{M}$
$\displaystyle+\left[\left(\bar{\bar{Z}}_{R2}+\bar{\bar{Z}}_{C2}+\bar{\bar{Z}}_{RL}+\bar{\bar{Z}}_{L2}-\bar{\bar{Z}}_{M}\right)^{-1}+\bar{\bar{Z}}_{M}^{-1}\right]^{-1}.$

The current can be calculated as the dot multiplication of the admittance
matrix (inverse of the impedance matrix) and the voltage vector:

$\mathbf{i}_{s}=\bar{\bar{Z}}_{\rm tot}^{-1}\cdot\mathbf{v}_{s}.$ (11)

Similarly, we can calculate the voltage and current in the load impedance,
$\mathbf{i}_{L}$ and $\mathbf{v}_{L}$:

$\mathbf{i}_{L}=(\bar{\bar{Z}}_{L2}-\bar{\bar{Z}}_{M}+\bar{\bar{Z}}_{R2}+\bar{\bar{Z}}_{C2}+\bar{\bar{Z}}_{RL})^{-1}\cdot\mathbf{v}_{M},$
(12)

where
$\mathbf{v}_{M}=\mathbf{v}_{s}-(\bar{\bar{Z}}_{s}+\bar{\bar{Z}}_{C1}+\bar{\bar{Z}}_{R1}+\bar{\bar{Z}}_{L1}-\bar{\bar{Z}}_{M})\cdot\mathbf{i}_{s}$.
Then, the voltage across the load can be found as
$\mathbf{v}_{L}=\bar{\bar{Z}}_{RL}^{-1}\cdot\mathbf{i}_{L}$. Knowing all the
currents and voltages in the circuit, we can calculate the efficiency by Eq.
(1)

##  III Efficiency analysis of time-modulated WPT systems

In this section, we analyse the effects of time-modulation of inductances on
the performance of WPT systems when the modulation is made at twice the
frequency of the power source. Here, we consider modulation of the self-
inductances of the Tx and Rx coils (i.e., $L_{1}(t)$ and $L_{2}(t)$) and of
the mutual inductance (i.e., $M(t)$) as expressed in the following relations:

$\displaystyle L_{\\{1,2\\}}(t)$ $\displaystyle=$ $\displaystyle
L_{0}\Big{[}1+m_{\rm L}{{}_{\\{1,2\\}}}\cos(\omega_{\rm
M}t+\phi_{\rm{L}}{{}_{\\{1,2\\}}})\Big{]}$ (13) $\displaystyle M(t)$
$\displaystyle=$ $\displaystyle M_{0}\Big{[}1+m_{\rm M}\cos(\omega_{\rm
M}t+\phi_{\rm{M}})\Big{]}.$ (14)

Here, $m_{L\\{1,2\\}}$ and $m_{\rm M}$ are the modulation depths of the self-
inductance and mutual inductance, respectively, and $\phi_{L\\{1,2\\}}$ and
$\phi_{\rm M}$ are the modulation phases of self and mutual inductances,
respectively. These parameters in principle can be controlled independently.
Here, for simplicity of the analysis, we consider the self and mutual
inductances to be modulated with the same modulation depth and phase, i.e.,
$m_{L\\{1,2\\}}=m_{\rm M}=m$ and $\phi_{L\\{1,2\\}}=\phi_{\rm M}=\phi$. For
this modulation case, we study possibilities of using such time-varying
circuits for improvement of the device performance. As an example, we assume
$Z_{\rm s}=2\leavevmode\nobreak\ \Omega$, $R_{L}=20\leavevmode\nobreak\
\Omega$, $R_{1}=R_{2}=0.1\leavevmode\nobreak\ \Omega$, $C_{1}=C_{2}=25.33$ nF,
$L_{0}=100\leavevmode\nobreak\ \mu$H, $\omega_{\rm s}=2\pi\times 100$ kHz,
$\omega_{\rm M}=2\omega_{\rm s}$. Figure 3 shows the efficiency and the power
received by the load with respect to the coupling coefficient $k$, for two
cases, with and without modulation. It can be seen that with modulation the
wireless transfer efficiency can be improved significantly. For small coupling
coefficients ($k\approx 0.01$), the efficiency is near zero when the coil is
not modulated. Strikingly, with modulation, the efficiency improves
significantly ($\eta>0.6$). It is also clear that the efficiency for the
modulated case is always larger than for the corresponding static system in a
large range of $k$. The increase in efficiency is accompanied with an increase
of the power transferred to the load. One can see that the load power is
always larger than that in the static system, i.e., $P_{L,m}/P_{L,s}>1$ for
all values of $k$ where $P_{L,m}$ and $P_{L,s}$ represent the received power
at the load for modulated and static cases, respectively. The results in Fig.
3(a) are obtained via optimization of $m$ and $\phi$. The optimized modulation
parameters for each value of $k$ are displayed in Fig. 3(b). It should be
noted that the improvement of efficiency under time modulation is not caused
by the enhanced mutual inductance within one half of the modulation period,
$M_{0}(1+m)$. This can be simply seen in Fig. 3(a), where the efficiency of a
static WPT system with $M=M_{0}(1+m)$ is plotted (see black dotted line).

Figure 3: (a) Power transfer efficiency and the output power at the load with
and without modulation, in terms of the coupling coefficient. The modulation
parameters for each value of $k$ are optimized independently. (b) The
optimized modulation parameters for each value of $k$.

Figure 4: (a) Optimized efficiencies, (b) modulation depths, and (c)
modulation phases for different values of the load resistance.

Next, we investigate the impact of individual component modulation on the
efficiency of the WPT system. We begin by examining the effect of modulating
$L_{1}$, which induces an effective negative resistance in the primary branch
at the optimal phase angle [28]. Basically, in this case the time-modulated
inductance acts as a parametric amplifier in the transmitter circuit. This
results in an increase of the current in the transmitter and receiver loops,
thereby boosting the output power. However, although the modulated $L_{1}$
configuration yields higher power, it does not necessarily enhance efficiency.
In fact, the modulation of only $L_{1}$ increases losses in proportion to the
power increase, thus, no efficiency improvement is achieved. On the other
hand, modulating $L_{2}$ emulates a negative resistance in the receiver loop,
resulting in an increase in the power delivered to the receiver without
incurring additional losses in the transmitter coil or in the source
resistance. This leads to a significant improvement of efficiency. Similarly,
the modulation of $M$ pumps energy directly to the receiver loop at the
optimal phase angle, which also contributes to the efficiency enhancement of
the WPT system. Obviously, there are inevitable additional losses in the
modulation circuits, but since these circuits do not carry significant power
(they only modulate the reactive components), these additional losses are
expected to be not significant in comparison to dissipation in the coils and
in the main power source.

In order to evaluate the effects of time modulations of separate inductors, we
fix the coupling coefficient to $k=0.1$ and examine the efficiency improvement
as a function of the load resistance (presented in Fig. 4). Here, we consider
three cases, modulating only $M$, $L_{2}$, and $M$ and $L_{2}$ simultaneously.
Figure 4(a) shows the efficiency in these scenarios. Note that the modulation
phase is optimized to achieve maximum efficiency at each load. We see that the
mutual inductance modulation (modulating only $M$) has a positive effect on
the efficiency improvement, for all values of the considered load resistance
($1\leavevmode\nobreak\ \Omega<R_{L}<1000\leavevmode\nobreak\ \Omega$).
Modulation of only $L_{2}$ has a much more positive effect. Therefore,
modulating $L_{2}$ and $M$ simultaneously (with the same modulation phase) is
potentially the optimal way for efficiency improvement. The modulation phases
and amplitudes for the three modulation scenarios are displayed in Fig. 4(b,c)
that are obtained from optimization. For all three cases, we confirm that the
received power at the load is always larger than in the static case.

##  IV Experimental study

###  IV.1 Experiment setup

It is important to note that practical implementation of any electrically-
tunable inductance is very challenging. There are some attempts for creation
of tunable inductors based on MEMS [30], magnetoelectric materials [31], and
non-Foster-based systems [32]. However, all these systems are still in the
very early phase of the research. To the best of our knowledge, there are no
works on tuning mutual inductance. Therefore, for the experimental
verification of our concept we used a simple mechanical system to modulate the
mutual inductance between two coils. Specifically, a wireless power transfer
demonstrator based on a low-frequency mechanical realization of a time-varying
mutual inductance has been designed (Fig. 5). We focus our efforts on
realization of time-modulated mutual inductance, because modulation of
ordinary inductance creates better understood parametric amplification, while
it appears that the effects of modulations of coupling have not been studied
before.

A schematic diagram of the system (see Fig. 5(a)) resembles a classical WPT
system that transfers energy from the source ($V_{\rm osc}$) to a load
($R_{L}$) via a resonant pair of inductively coupled coils, forming a
transformer. However, instead of having a fixed mutual inductance (and fixed
transformation ratio), this transformer is a time-dependent device. It is
based on a time-varying mutual inductance controlled by a modulation signal
$V_{x}$ (a pump), as is illustrated in Fig. 5(b).

Figure 5: Experimental WPT demonstrator based on a mechanically time-varying
mutual inductance. (a) Schematic diagram. (b) Simplified construction. (c)
Modulation function of $L_{1}$, $L_{2}$, and $M$ in one temporal period.
$T_{\rm M}=2\pi/\omega_{\rm M}$ is the modulation period.

Figure 6: Photo of the experimental setup.

The transformer comprises two commercial E52 ferrite cores with a mutual air
gap. The primary and secondary coils have 62 turns of 0.6 mm diameter copper
wire, wound in a bifilar fashion. The coils were wound around both central
cores, one on top of the other. One side of the ferrite core is fixed to a
solid base, while the second one is mounted onto an electromechanical
transducer. The mechanical construction of the ‘moving joint’ was very
challenging since even very small mechanical vibrations (in the order of a
tenth of a millimetre) can affect the length of the air gap and, therefore,
the mutual inductance. For this reason, a precise two-axes mechanical
positioner was used for accurate adjustment of desired core locations. The
transducer is driven by a pump signal (the first output of a two-channel
arbitrary waveform generator, AWG) amplified with the help of an in-house two-
channel class D power amplifier. Thus, mechanical modulation of the mutual
inductance is achieved by movement of one of the ferrite cores, i.e., by
varying the length of the air gap. The second output of the arbitrary waveform
generator was connected to the second channel of the class D power amplifier,
feeding the primary coil of a transformer while the secondary coil was loaded
by a resistive load. Thanks to high-efficient D-class operation, no
significant thermal problems were noted even with very high output current
levels (up to 3 A), needed for excitation of the mechanical transducer.

The circuit parameters in the experiment are as following:
$R_{1}=R_{2}=0.3\leavevmode\nobreak\ \Omega$, $C_{1}=C_{2}=88$ mF,
$R_{L}=4.7\leavevmode\nobreak\ \Omega$, $V_{\rm osc}=5$ V with the source
resistance $Z_{\rm s}=5.6\leavevmode\nobreak\ \Omega$. The modulation (pump)
frequency is limited by the mechanical self-resonance of the system, which is
about 110 Hz. In order to avoid parasitic effects caused by the mechanical
self-resonance, the modulation frequency should be lower than the self-
resonance frequency. Therefore, we chose the modulation frequency of
$\omega_{\rm M}=2\pi\times 90$ Hz and the source frequency of $\omega_{\rm
s}=2\pi\times 45$ Hz. The modulation profiles for $L_{1}(t)$, $L_{2}(t)$, and
$M(t)$ are shown in Fig. 5(c). Furthermore, we paid special attention to
achieving the required phase shift between the input and pump signals. The
problem is that the phase shift set on the AWG does not match the actual phase
shift between the two signals, due to different delays in the two channels of
the power amplifier and the different lengths of the connecting cables.
Therefore, we set a zero phase at the AWG and checked the phase difference
between the signals at the transducer and at the primary winding of a
transformer with an oscilloscope. Then we introduced an additional phase shift
required to compensate for the delay difference. This correction is very
effective for the referent phase shift (zero degree), while its accuracy
decreases for other phase shifts. Another practical problem is the unknown
output impedance of the power amplifier (in the theoretical part this
impedance was assumed to be zero). We do not know the exact value of this
impedance, but a “rule of thumb” estimate gives a value of about 0.1 $\Omega$
to 0.3 $\Omega$. Therefore, we calculated the associated voltage attenuation
ratios and compensated for them by making appropriate adjustments to the gains
of the oscilloscope channels. After all these initial calibration procedures
were completed, we began the measurements that will be explained in the next
section.

Figure 7: (a)(b)(c) The theoretical (orange line) and measured (blue line)
voltages at the primary coil $V_{1}$ and (d)(e)(f) the secondary coil $V_{2}$
for the modulation phases of $\phi=9^{\circ}$, $\phi=54^{\circ}$, and
$\phi=99^{\circ}$.

###  IV.2 Results

The theoretical and measured waveforms of $V_{1}$ and $V_{2}$ are compared in
Fig. 7. We consider three scenarios with the modulation phases
$\phi=9^{\circ}$, $\phi=54^{\circ}$, and $\phi=99^{\circ}$. It can be seen
that in all three scenarios the voltage profile in the time domain is not
sinusoidal, which is a clear evidence that higher-order harmonics are
generated. The measured signals are quite similar to those obtained from the
theory. There are slight differences that are more pronounced at larger phase
delays. This is a direct consequence of the limitations of the experimental
system described in the previous section. Nevertheless, the obtained results
clearly support the conclusions from the developed theory. This is also
evident from the table showing the measured power received by the load with
and without modulation. It can be seen that the efficiency is significantly
higher for some modulation phases than for the static system, for example,
$\phi=54^{\circ}$ and $\phi=99^{\circ}$ in our experiment. Of course, like in
any parametric system, this time-varying system is phase dependent. One must
synchronize the phases of the signal and modulation in order to obtain an
improved efficiency.

WPT Efficiency
---
Modulation phase | with modulation | without modulation
$\phi=9^{\circ}$ | 8.13% | 9.28%
$\phi=54^{\circ}$ | 19.8%
$\phi=99^{\circ}$ | 27.73%
Table 1: Measured efficiency with and without modulation for three different
modulation phases.

##  V Conclusion

In summary, this paper introduces an inductive wireless power transfer system
with coupling inductance varying periodically over time. The study finds that
by varying the inductance of the coupled coils and the mutual inductance, the
WPT system’s efficiency can be increased significantly. At the same time, the
power received by the load is also improved as compared to the reference
static WPT system. This technique overcomes the fundamental tradeoff between
efficiency and delivered power, improving the transferred power and transfer
efficiency simultaneously. To validate this conclusion, a conceptual
experiment was conducted, demonstrating that modulation of the coil inductance
and mutual inductance can lead to a significant improvement in power transfer
efficiency. The experimental results confirm the theoretical conclusions drawn
from the study.

## References

  * [1] W. Brown, “The history of power transmission by radio waves,” _IEEE Transactions on Microwave Theory and Techniques_ , vol. 32, no. 9, pp. 1230–1242, 1984.
  * [2] A. Kurs, A. Karalis, R. Moffatt, J. D. Joannopoulos, P. Fisher, and M. Soljac̆ić, “Wireless power transfer via strongly coupled magnetic resonances,” _Science_ , vol. 317, no. 5834, pp. 83–86, 2007.
  * [3] S. Li and C. C. Mi, “Wireless power transfer for electric vehicle applications,” _IEEE Journal of Emerging and Selected Topics in Power Electronics_ , vol. 3, no. 1, pp. 4–17, 2015.
  * [4] M. Song, P. Jayathurathnage, E. Zanganeh, M. Krasikova, P. Smirnov, P. Belov, P. Kapitanova, C. Simovski, S. Tretyakov, and A. Krasnok, “Wireless power transfer based on novel physical concepts,” _Nature Electronics_ , vol. 4, pp. 707–716, 2021.
  * [5] S. Y. R. Hui, W. Zhong, and C. K. Lee, “A critical review of recent progress in mid-range wireless power transfer,” _IEEE Transactions on Power Electronics_ , vol. 29, no. 9, pp. 4500–4511, 2013.
  * [6] A. P. Sample, D. T. Meyer, and J. R. Smith, “Analysis, experimental results, and range adaptation of magnetically coupled resonators for wireless power transfer,” _IEEE Transactions on Industrial Electronics_ , vol. 58, no. 2, pp. 544–554, 2011.
  * [7] W. Zhang and C. C. Mi, “Compensation topologies of high-power wireless power transfer systems,” _IEEE Transactions on Vehicular Technology_ , vol. 65, no. 6, pp. 4768–4778, 2016.
  * [8] Q. Li and Y. C. Liang, “An inductive power transfer system with a high-Q resonant tank for mobile device charging,” _IEEE Transactions on Power Electronics_ , vol. 30, no. 11, pp. 6203–6212, 2015.
  * [9] T. P. Duong and J.-W. Lee, “Experimental results of high-efficiency resonant coupling wireless power transfer using a variable coupling method,” _IEEE Microwave and Wireless Components Letters_ , vol. 21, no. 8, pp. 442–444, 2011.
  * [10] J. P. K. Sampath, A. Alphones, D. M. Vilathgamuwa, A. Ong, and X. B. Nguyen, “Coil enhancements for high efficiency wireless power transfer applications,” in _IECON 2014 - 40th Annual Conference of the IEEE Industrial Electronics Society_ , pp. 2978–2983, 2014.
  * [11] S. Cheon, Y.-H. Kim, S.-Y. Kang, M. L. Lee, J.-M. Lee, and T. Zyung, “Circuit-model-based analysis of a wireless energy-transfer system via coupled magnetic resonances,” _IEEE Transactions on Industrial Electronics_ , vol. 58, no. 7, pp. 2906–2914, 2011.
  * [12] S. R. Hui, “Magnetic resonance for wireless power transfer [a look back],” _IEEE Power Electronics Magazine_ , vol. 3, no. 1, pp. 14–31, 2016.
  * [13] Q. Deng, J. Liu, D. Czarkowski, M. K. Kazimierczuk, M. Bojarski, H. Zhou, and W. Hu, “Frequency-dependent resistance of litz-wire square solenoid coils and quality factor optimization for wireless power transfer,” _IEEE Transactions on Industrial Electronics_ , vol. 63, no. 5, pp. 2825–2837, 2016\.
  * [14] A. L. F. Stein, P. A. Kyaw, and C. R. Sullivan, “High-Q self-resonant structure for wireless power transfer,” in _2017 IEEE Applied Power Electronics Conference and Exposition (APEC)_ , pp. 3723–3729, 2017.
  * [15] Y. Zhang and Z. Zhao, “Frequency splitting analysis of two-coil resonant wireless power transfer,” _IEEE Antennas and Wireless Propagation Letters_ , vol. 13, pp. 400–402, 2014.
  * [16] K. Agarwal, R. Jegadeesan, Y.-X. Guo, and N. V. Thakor, “Wireless power transfer strategies for implantable bioelectronics,” _IEEE reviews in biomedical engineering_ , vol. 10, pp. 136–161, 2017.
  * [17] A. Zakerian, S. Vaez-Zadeh, and A. Babaki, “A dynamic WPT system with high efficiency and high power factor for electric vehicles,” _IEEE Transactions on Power Electronics_ , vol. 35, no. 7, pp. 6732–6740, 2019.
  * [18] D.-W. Seo, J.-H. Lee, and H.-S. Lee, “Optimal coupling to achieve maximum output power in a WPT system,” _IEEE Transactions on Power Electronics_ , vol. 31, no. 6, pp. 3994–3998, 2015.
  * [19] F. Mastri, A. Costanzo, and M. Mongiardo, “Coupling-independent wireless power transfer,” _IEEE Microwave and Wireless Components Letters_ , vol. 26, no. 3, pp. 222–224, 2016.
  * [20] E. Galiffi, R. Tirole, S. Yin, H. Li, S. Vezzoli, P. A. Huidobro, M. G. Silveirinha, R. Sapienza, A. Alù, and J. Pendry, “Photonics of time-varying media,” _Advanced Photonics_ , vol. 4, no. 1, pp. 014 002–014 002, 2022.
  * [21] G. Ptitcyn, M. S. Mirmoosa, A. Sotoodehfar, and S. A. Tretyakov, “A tutorial on the basics of time-varying electromagnetic systems and circuits: Historic overview and basic concepts of time-modulation.” _IEEE Antennas and Propagation Magazine_ , 2023.
  * [22] Z. Hayran and F. Monticone, “Using time-varying systems to challenge fundamental limitations in electromagnetics: Overview and summary of applications.” _IEEE Antennas and Propagation Magazine_ , 2023.
  * [23] D. L. Sounas and A. Alù, “Non-reciprocal photonics based on time modulation,” _Nature Photonics_ , vol. 11, no. 12, pp. 774–783, 2017.
  * [24] Z. Wu and A. Grbic, “Serrodyne frequency translation using time-modulated metasurfaces,” _IEEE Transactions on Antennas and Propagation_ , vol. 68, no. 3, pp. 1599–1606, 2019.
  * [25] X. Wang, M. S. Mirmoosa, V. S. Asadchy, C. Rockstuhl, S. Fan, and S. A. Tretyakov, “Metasurface-based realization of photonic time crystals,” _Science Advances_ , vol. 9, no. 14, p. eadg7541, 2023.
  * [26] H. Li, A. Mekawy, and A. Alù, “Beyond Chu’s limit with floquet impedance matching,” _Phys. Rev. Lett._ , vol. 123, p. 164102, Oct. 2019.
  * [27] X. Yang, E. Wen, and D. F. Sievenpiper, “Broadband time-modulated absorber beyond the Bode-Fano limit by energy trapping,” _arXiv preprint arXiv:2110.12843_ , 2021.
  * [28] P. Jayathurathnage, F. Liu, M. S. Mirmoosa, X. Wang, R. Fleury, and S. A. Tretyakov, “Time-varying components for enhancing wireless transfer of power and information,” _Physical Review Applied_ , vol. 16, no. 1, p. 014017, 2021\.
  * [29] J. Bird, _Electrical Circuit Theory and Technology_ , 6th ed., ch. 46, p. 736\. New York: Routledge, 2017.
  * [30] I. El Gmati, P. Calmon, A. Boukabache, P. Pons, H. Boussetta, M. Kallala, and K. Besbes, “Liquid RF MEMS variable inductor,” _Procedia Engineering_ , vol. 5, pp. 1380–1383, 2010.
  * [31] J. Lou, D. Reed, M. Liu, and N. Sun, “Electrostatically tunable magnetoelectric inductors with large inductance tunability,” _Applied Physics Letters_ , vol. 94, no. 11, 2009.
  * [32] E. Vorobev, V. Turgaliev, D. Kholodnyak, and N. Ivanov, “Active tunable inductor using non-Foster element,” in _2017 IEEE Conference of Russian Young Researchers in Electrical and Electronic Engineering (EIConRus)_ , pp. 346–349. IEEE, 2017.

*

##  Appendix A Validity of the T-circuit interpretation for time-varying
coupled coils

It is known that two inductively coupled coils can be considered equivalent to
a T-circuit with inductances $L_{1}-M$, $L_{2}-M$, and $M$. This equivalence
is usually derived for static components, meaning that $L_{1}$, $L_{2}$, and
$M$ are time-invariant. Next, we demonstrate that the T-circuit interpretation
of mutual inductance is also valid when the components are varying as
arbitrary functions of time.

In the time domain, the voltage at the two sides of a coil (see Fig. 2(top))
can be expressed as

$v_{1}(t)=\frac{d[L_{1}(t)i_{1}(t)]}{dt}+\frac{d[M(t)i_{2}(t)]}{dt},$ (15)
$v_{2}(t)=\frac{d[L_{2}(t)i_{2}(t)]}{dt}+\frac{d[M(t)i_{1}(t)]}{dt}.$ (16)

These two equations can be rearranged as

$v_{1}(t)=\frac{d\left\\{[L_{1}(t)-M(t)]i_{1}(t)\right\\}}{dt}+\frac{d\left\\{M(t)[i_{1}(t)+i_{2}(t)]\right\\}}{dt},$
(17)
$v_{2}(t)=\frac{d\left\\{[L_{2}(t)-M(t)]i_{2}(t)\right\\}}{dt}+\frac{d\left\\{M(t)[i_{1}(t)+i_{2}(t)]\right\\}}{dt}.$
(18)

One can see that the above two equations are exactly the voltage-current
relation of a T-circuit (see Fig. 2(bottom)). Therefore, time-varying
magnetically coupled coils can still be considered equivalent to a T-circuit.
This property brings convenience for the analysis of time-varying coils.