#  Effect of collective spin excitation on electronic transport in topological
spin texture

Kohei Hattori<EMAIL_ADDRESS>Department of Applied
Physics, The University of Tokyo, Bunkyo, Tokyo 113-8656, Japan  Hikaru
Watanabe  0000-0001-7329-9638<EMAIL_ADDRESS>Research
Center for Advanced Science and Technology, The University of Tokyo, Meguro,
Tokyo 153-8904, Japan  Junta Iguchi  Department of Applied Physics, The
University of Tokyo, Bunkyo, Tokyo 113-8656, Japan  Takuya Nomoto
0000-0002-4333-6773 Research Center for Advanced Science and Technology, The
University of Tokyo, Meguro, Tokyo 153-8904, Japan  Department of Physics,
Tokyo Metropolitan University, Hachioji, Tokyo 192-0397, Japan  Ryotaro Arita
0000-0001-5725-072X Research Center for Advanced Science and Technology, The
University of Tokyo, Meguro, Tokyo 153-8904, Japan  Center for Emergent Matter
Science, RIKEN, Wako, Saitama 351-0198, Japan

###### Abstract

We develop an efficient real-time simulation method for the spin-charge
coupled system in the velocity gauge. This method enables us to compute the
real-time simulation for the two-dimensional system with the complex spin
texture. We focus on the effect of the collective excitation of the localized
spins on the electronic transport properties of the non-trivial topological
state in real space. To investigate this effect, we calculate the linear
optical conductivity by calculating the real-time evolution of the Kondo
lattice model on the triangular lattice, which hosts an all-in/all-out
magnetic structure. In the linear conductivity spectra, we observe multiple
peaks below the bandgap regime, attributed to the resonant contributions of
collective modes similar to the skyrmionic system, alongside broadband
modifications resulting from off-resonant spin dynamics. The result shows that
the collective excitation, similar to the skyrmionic system, influences the
optical response of the electron systems based on symmetry analysis. We
elucidate the interference between the contributions from the different spin
excitations to the optical conductivity in the multiple spin texture, pointing
out the mode-dependent electrical activity. We show the complex interplay
between the complex spin texture and the itinerant electrons in the two-
dimensional spin-charge coupled system.

##  I introduction

Recent research intensively elucidated the effect of collective dynamics of
the spontaneous order, such as magnetic order and excitonic order, on the
optical responses of electronic systems [1, 2, 3, 4, 5, 6, 7]. For instance,
theoretical studies simultaneously calculate the real-time evolution of the
electron system and the order parameter [4, 5, 6, 7]. More specifically, the
study of magnetic order [4] shows that the collective excitation of the
localized spins modulates the linear optical conductivity and photocurrent
response in an antiferromagnetic chain. The results identified the features
which are absent within the independent particle approximation. It is,
however, not straightforward to apply the methodology to more complex systems,
such as those in a two-dimensional and with a complex spin structure, due to
their large computational cost. The realization of efficient real-time
simulation applicable to diverse cases is highly desirable since it is
anticipated that one can explore intriguing spin-charge coupled dynamics, for
example, by considering the magnetic order whose real-space topological
texture may significantly influence electronic properties [8, 9, 10, 11].

The non-trivial topological state in real space appears in non-coplanar
magnetic structures such as frustrated magnets and magnetic skyrmionic
systems, where the complex spin texture formed by spin-orbit coupling [12] and
itinerant-electron-mediated interactions [13, 14]. Through the exchange
interaction with localized spin systems having such nontrivial spin textures,
the electron system feels the fictitious magnetic field and exhibits anomalous
transverse transport. Theoretical studies [8, 9, 15, 16, 17, 18, 19] have
shown that the anomalous transverse conduction occurs in the Kondo lattice
model with nontrivial spin structures. Importantly, such non-coplanar spin
textures, including the all-in/all-out (AIAO) structure [20, 21, 22] and
skyrmion crystal [13, 23, 14], allows for collective modes richer than those
in conventional ferro- and antiferromagnets. Indeed, theoretical [24, 25, 26]
and experimental studies [27, 28, 29] identified collective modes of the
localized spins in the skyrmionic system. These collective modes include
breathing modes excited by out-of-plane magnetic fields and clockwise and
counterclockwise rotatory modes excited by in-plane magnetic fields [24, 25].

In this study, we elaborate on the effect of the collective excitation of the
AIAO spin texture on the anomalous transverse transport of the electron
system. We simulate the real-time evolution of the Kondo lattice model on the
triangular lattice hosting the AIAO spin texture comprised of four sites in
the unit cell. The dynamical property is clarified by computing the real-time
evolution of the electrons and localized spins simultaneously. The
computational cost of the two-dimensional spin-charge coupled system with the
AIAO spin texture is expensive because of the dimensionality and complex spin
structure. We resolve the difficulty by developing an efficient real-time
simulation based on the velocity-gauge formulation and interpolation
technique. The developed calculation scheme allows us to demonstrate the
intriguing interplay between the complex spin texture and itinerant electrons
in a two-dimensional system.

The real-time simulation unambiguously shows that the optical conductivity
spectrum of the electron system exhibits multiple peaks originating from rich
collective dynamics of localized spins. To classify which collective mode
affects linear optical conductivity, we analyze the spin excitations by the
symmetry-adapted bases. The symmetry analysis indicates the similarity between
the collective modes of the AIAO spin system and those of the skyrmionic
system. We find that the rotatory modes of the localized spins affect the
optical responses of the electronic systems in two-dimensional systems.
Furthermore, the interference effect of the collective modes, which is
characteristic of the system hosting the complex spin structure, results in
the in-gap spectrum of the optical Hall conductivity whose intensity
significantly differs between the modes. These findings may lead to an
understanding of the interplay between the complex spin texture and itinerant
electrons.

The paper is organized as follows. In Sec. II, we explain the details of the
model and the computational scheme for the real-time evolution of the spin-
charge coupled system. Section III is for the symmetry classification of the
spin excitations and its coupling to the external stimuli. Then, we elucidate
the physical responses and corroborate the effects of the spin-charge coupled
dynamics in Sec. IV. We draw the conclusion in Sec. V.

##  II method

###  II.1 Model

Figure 1: (a) Kondo lattice model with AIAO order on the triangular lattice
with four-sublattice unit cell. (b) Spin moments at each sublattice. The spin
texture is the AIAO state including canted component along the $z$-axis.

We work on the AIAO state in the Kondo lattice model on the triangular lattice
(FIG. 1). The tight-binding Hamiltonian is expressed as

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=-\sum_{\langle
ij\rangle\sigma\sigma^{\prime}}t_{h}\hat{c}_{i\sigma}^{{\dagger}}\hat{c}_{j\sigma^{\prime}}-J\sum_{i\sigma\sigma^{\prime}}\hat{c}_{i\sigma}^{{\dagger}}\bm{\mathrm{S}}_{i}\cdot\bm{\sigma}_{\sigma\sigma^{\prime}}\hat{c}_{i\sigma^{\prime}}$
$\displaystyle\quad-i\lambda\sum_{\langle
ij\rangle\sigma\sigma^{\prime}}\hat{c}_{i\sigma}^{{\dagger}}\bm{e}_{ij}\cdot\bm{\sigma}_{\sigma\sigma^{\prime}}\hat{c}_{j\sigma^{\prime}}-K_{z}\sum_{i}({\mathrm{S}}_{i}^{z})^{2}.$
(1)

$\hat{c}_{i\sigma}^{\dagger}$ ($\hat{c}_{i\sigma}$) is the creation
(annihilation) operator of the electron on the site $i$ with the spin index
$\sigma$ ($\sigma=\uparrow,\downarrow$). The first term of the Hamiltonian
represents the hopping of the electrons between the nearest neighbor site
$\langle ij\rangle$ with the amplitude $t_{h}$. The second term represents the
exchange interaction $J$ between the spin moment $\bm{\sigma}$ of the
itinerant electron and the localized spin moment $\bm{\mathrm{S}}$. The third
term represents the antisymmetric spin-orbit coupling $\lambda$ (ASOC), which
acts as an effective magnetic field parallel to the vector $\bm{e}_{ij}$
between the nearest neighboring sites $\langle ij\rangle$. The ASOC breaks the
inversion symmetry and thereby allows the system to exhibit the cross
correlation between the spin and charge degrees of freedom. The last term
represents the easy anisotropy of the localized spin system $K_{z}$.

In the following calculations we use the parameters
$t_{h}=1,~{}J=3.0,~{}\lambda=0.3,~{}K_{z}=0.1,~{}\alpha_{G}=0.01,~{}\gamma=0.01$,
unless explicitly mentioned. The filling number $n$ of electrons is set to
$0.25$ per unit cell to stabilize the AIAO state in equilibrium [FIG. 1 (b)]
[15, 30], where the system shows insulating behavior. Owing to the anisotropy
of the spin Hamiltonian, a finite magnetization appears along the $z$-axis. We
also set lattice constant $a=1$ and elementary charge $e=1$.

We transform the Hamiltonian $\hat{\mathcal{H}}$ in the real-space
representation into the momentum representation through the Fourier
transformation

$\displaystyle\hat{c}^{{\dagger}}_{a\sigma}(\bm{k})=\sum_{\bm{R}}e^{i\bm{k}\cdot(\bm{r}_{a}+\bm{R})}\hat{c}^{{\dagger}}_{a^{\prime}\sigma},$
(2)

where the translation vectors $\bm{R}$ are defined as
$\bm{R}=n_{1}(2,0,0)+n_{2}(1,\sqrt{3},0)$ $(n_{1},n_{2}\in\mathbb{Z})$ and the
position $\bm{r}_{a^{\prime}}$ of the lattice site $a^{\prime}$ is
$\bm{r}_{a^{\prime}}=\bm{r}_{a}+\bm{R}$. Since the AIAO state has the unit
cell four times larger than that in the paramagnetic state, the index $a$ for
the sublattice in $\hat{c}^{{\dagger}}_{a\sigma}(\bm{k})$ is taken as
$a=1,2,3,4$ [FIG. 1 (a)]. The Hamiltonian $\hat{\mathcal{H}}$ in the momentum
representation reads as

$\displaystyle\hat{\mathcal{H}}=\sum_{\bm{k}}\sum_{\sigma\sigma^{\prime}}\sum_{ab}\quantity(\bm{\mathcal{H}}_{0}(\bm{k}))^{\sigma\sigma^{\prime}}_{ab}\hat{c}_{a\sigma}^{\dagger}(\bm{k})\hat{c}_{b\sigma^{\prime}}(\bm{k}),$
(3)

where $\bm{\mathcal{H}}_{0}(\bm{k})$ is the electronic Hamiltonian at each
$k$-point.

When the electric field is incorporated under the velocity gauge, the creation
(annihilation) operator changes as $\hat{c}_{i\sigma}^{{\dagger}}\rightarrow
e^{i\bm{A}(t)\cdot\bm{r}}\hat{c}_{i\sigma}^{{\dagger}}\
(\hat{c}_{i\sigma}\rightarrow e^{-i\bm{A}(t)\cdot\bm{r}}\hat{c}_{i\sigma})$,
where $\bm{A}(t)$ is the vector potential of the electromagnetic field. In the
velocity gauge, a time-dependent Hamiltonian is

$\displaystyle\hat{\mathcal{H}}(t)=\sum_{\bm{k}}\sum_{\sigma\sigma^{\prime}}\sum_{ab}\quantity(\bm{\mathcal{H}}(\bm{k},t))^{\sigma\sigma^{\prime}}_{ab}\hat{c}_{a\sigma}^{\dagger}(\bm{k})\hat{c}_{b\sigma^{\prime}}(\bm{k}),$
(4)

where the time-dependent electronic Hamiltonian $\bm{\mathcal{H}}(\bm{k},t)$
at each $k$-point is expressed as

$\displaystyle\bm{\mathcal{H}}(\bm{k},t)=\bm{\mathcal{H}}_{0}(\bm{k}-\bm{A}(t)),$
(5)

while the creation and annihilation operators are coincident with the original
operators in Eq. (3) without the dependence on $t$. The light-matter coupling
is given by the time-dependent vector potential $\bm{A}(t)$ in the velocity
gauge, which indicates the photo-electric field
$\bm{E}(t)=-\partial_{t}\bm{A}(t)$. Since the light-matter coupling is taken
into account with the electric-dipole approximation, the vector potential is
spatially uniform and thus does not break the translation symmetry. In the
expression Eq. (4), we assume that the Wannier state of an electron is well-
localized at a given site.

###  II.2 Calculation Scheme

In this study, we calculate the real-time evolution of the spin-charge coupled
system. To compute the real-time simulation, we simultaneously solve the von
Neumann equation and the Landau-Lifshitz-Gilbert (LLG) equation.

First, the time evolution of electrons can be described by the single-particle
density matrix (SPDM)
$\rho^{\sigma\sigma^{\prime}}_{ab}(\bm{k},t)=\langle\hat{c}^{{\dagger}}_{b\sigma^{\prime}}(\bm{k})\hat{c}_{a\sigma}(\bm{k})\rangle_{t}\equiv\mathrm{Tr}[\hat{c}^{{\dagger}}_{b\sigma^{\prime}}(\bm{k})\hat{c}_{a\sigma}(\bm{k})\hat{\rho}(t)],$
where $\hat{\rho}(t)$ is the density matrix of the system at time $t$. The
SPDM satisfies the following equation called the von Neumann equation [31],

$\displaystyle\frac{\partial\bm{\rho}(\bm{k},t)}{\partial
t}=-i[\bm{\mathcal{H}}(\bm{k},t),\bm{\rho}(\bm{k},t)]-\gamma(\bm{\rho}(\bm{k},t)-\bm{\rho}_{\mathrm{eq}}(\bm{k},t)).$
(6)

The vector potential $\bm{A}(t)$ for the light field does not break the
translation symmetry, which keeps the von Neumann equation block-diagonal with
respect to the momentum $\bm{k}$. The $\bm{k}$-local property, which is
compatible with the parallel computation, is advantageous for efficient
simulation of the real-time evolution of spin-charge coupled dynamics. On the
other hand, if one works on the length gauge under which the light field is
expressed by the scalar potential as Hamiltonian $-\bm{E}(t)\cdot\bm{r}$, the
translation-symmetry breaking leads to $k$-derivative of the SPDM
$-\bm{E}\cdot\frac{\partial\bm{\rho}(\bm{k},t)}{\partial\bm{k}}$ in the von
Neumann equation. This $k$-derivative term is unfavorable for parallel
computation of $k$-point grids.

Second, the time evolution of the localized spin system is governed by the LLG
equation

$\displaystyle\frac{d\bm{\mathrm{S}}_{a}}{dt}$
$\displaystyle=\frac{1}{1+\alpha_{\mathrm{G}}^{2}}\left(\bm{\mathrm{H}}_{a}^{\mathrm{eff}}\times\bm{\mathrm{S}}_{a}+\alpha_{\mathrm{G}}\bm{\mathrm{S}}_{a}\times(\bm{\mathrm{S}}_{a}\times\bm{\mathrm{H}}^{\mathrm{eff}}_{a})\right),$
(7) $\displaystyle\bm{\mathrm{H}}^{\mathrm{eff}}_{a}$
$\displaystyle=-J\langle\bm{\sigma}_{a}\rangle-2K_{z}(0,0,\mathrm{S}_{a}^{z}).$
(8)

Here, $\bm{\mathrm{S}}_{a}$ represents the localized spin moment at sublattice
site $a$, and $\bm{\mathrm{H}}_{a}^{\mathrm{eff}}$ represents the effective
magnetic field coupled to $\bm{\mathrm{S}}_{a}$. In the LLG equation Eq. (7),
$\left\langle{{\bf\it\sigma}_{a}}\right\rangle$ is the sublattice-dependent
spin density of electrons calculated with SPDM [see Eq. (13)]. Furthermore, we
account for a relaxation effect to obtain a physically reasonable response to
light. Though relaxation stems from electron-electron correlations, electron-
phonon interactions, and impurity scattering, we treat it phenomenologically
by using the relaxation time approximation in the von Neumann equation as
$\gamma({\bf\it\rho}(\bm{k},t)-{\bf\it\rho}_{\mathrm{eq}}(\bm{k},t))$ in Eq.
(6), and the Gilbert damping $\alpha_{G}$ in Eq. (7).
${\bf\it\rho}_{\text{eq}}(\bm{k},t)$ is the SPDM in equilibrium at the
temperature $T=0$. $\bm{\rho}_{\mathrm{eq}}(\bm{k},t)$ represents an
equilibrium property but shows the dependence on time $t$ when the light-field
is incorporated under the velocity gauge [32, 33] (discussed below). We solve
the coupled equations Eq. (6) and Eq. (7) by the fourth-order Runge-Kutta
method.

In the velocity gauge, the time-dependent Hamiltonian indicates the
modification of the velocity operator called the diamagnetic correction. The
correction effectively gives rise to the shift in momentum $\bm{k}$, and thus
the momentum-resolved SPDM in equilibrium differs between those at initial
time $t=0$ by the momentum shift. The momentum shift explicitly reads as

$\displaystyle\bm{\rho}_{\mathrm{eq}}(\bm{k},t)=\bm{\rho}^{0}_{\mathrm{eq}}(\bm{k}-\bm{A}(t)),$
(9)

where $\bm{\rho}^{0}_{\mathrm{eq}}(\bm{k})$ is the SPDM for the initial
Hamiltonian $\bm{\mathcal{H}}_{0}(\bm{k})$. The SPDM for the equilibrium
state, $\bm{\tilde{\rho}}_{\mathrm{eq}}(\bm{k},t)$ at time $t$, is calculated
in the Bloch basis as follows

$\displaystyle[\bm{\tilde{\rho}}_{\mathrm{eq}}(\bm{k},t)]_{mn}=\delta_{mn}\Theta(\mu-\epsilon_{n}(\bm{k},t)).$
(10)

We can compute the SPDM for the equilibrium state in the orbital basis by
using the unitary transformation

$\displaystyle\bm{\rho}_{\mathrm{eq}}(\bm{k},t)=\bm{U}(\bm{k},t)\bm{\tilde{\rho}}_{\mathrm{eq}}(\bm{k},t)\bm{U}^{{\dagger}}(\bm{k},t),$
(11)

where the unitary matrix $\bm{U}(\bm{k},t)$ diagonalizes the Hamiltonian
$\bm{\mathcal{H}}(\bm{k},t)$ as
$\bm{U}^{{\dagger}}(\bm{k},t)\bm{\mathcal{H}}(\bm{k},t)\bm{U}(\bm{k},t)=\bm{\mathcal{E}}(\bm{k},t)$
and $[\bm{\mathcal{E}}(\bm{k},t)]_{mn}=\delta_{mn}\epsilon_{n}(\bm{k},t)$. The
Hamiltonian changes as
$\bm{\mathcal{H}}(\bm{k},t)=\bm{\mathcal{H}}_{0}(\bm{k}-\bm{A}(t))$ due to the
electric field. Therefore, we need to calculate the SPDM
$\bm{\rho}_{\mathrm{eq}}(\bm{k},t)$ at every time step by diagonalization of
the Hamiltonian $\bm{\mathcal{H}}(\bm{k},t)$.

To avoid this large computational cost, we approximate the time-dependent SPDM
in equilibrium. We calculate the SPDM $\bm{\rho}_{\mathrm{eq}}(\bm{k},t)$ by
applying cubic interpolation of the momentum grid. We initially calculate the
SPDM $\bm{\rho}^{0}_{\mathrm{eq}}(\bm{k})$ at $t=0$ and obtain the
interpolated function $\bm{\rho}^{\mathrm{int}}_{\mathrm{eq}}(\bm{k})$ by
performing the cubic interpolation of $\bm{\rho}^{0}_{\mathrm{eq}}(\bm{k})$
with respect to $\bm{k}$. Then we can approximate the SPDM
$\bm{\rho}_{\mathrm{eq}}(\bm{k},t)$ at the time $t$ as

$\displaystyle\bm{\rho}_{\mathrm{eq}}(\bm{k},t)\simeq\bm{\rho}^{\mathrm{int}}_{\mathrm{eq}}(\bm{k}-\bm{A}(t)).$
(12)

By this interpolation technique, we can calculate
$\bm{\rho}_{\mathrm{eq}}(\bm{k},t)$ easily without diagonalization of the
Hamiltonian $\bm{\mathcal{H}}(\bm{k},t)$ at each time step.

Next, we explain a calculation method for physical responses. At each time
step, we evaluate the physical quantity $\langle\hat{\mathcal{O}}(t)\rangle$
as

$\displaystyle\langle\hat{\mathcal{O}}(t)\rangle=\frac{1}{N}\sum_{\bm{k}}\mathrm{Tr}[\bm{\rho}(\bm{k},t)\bm{\mathcal{O}}(\bm{k},t)],$
(13)

by using the SPDM. The current density operator is written by

$\displaystyle\hat{\bm{j}}(t)$
$\displaystyle=\frac{1}{n_{\mathrm{sub}}}\sum_{\bm{k}}\sum_{ab}\sum_{\sigma\sigma^{\prime}}\partialderivative{\quantity({\bf\it\mathcal{H}}(\bm{k},t))_{ab}^{\sigma\sigma^{\prime}}}{\bm{k}}\hat{c}_{a\sigma}^{\dagger}(\bm{k})\hat{c}_{b\sigma^{\prime}}(\bm{k}),$
$\displaystyle\equiv\frac{1}{n_{\mathrm{sub}}}\sum_{\bm{k}}\sum_{ab}\sum_{\sigma\sigma^{\prime}}\quantity({\bf\it
j}(\bm{k},t))_{ab}^{\sigma\sigma^{\prime}}\hat{c}_{a\sigma}^{\dagger}(\bm{k})\hat{c}_{b\sigma^{\prime}}(\bm{k}),$
(14)

where $n_{\mathrm{sub}}$ represents the number of sublattices in a unit cell.
Then, the linear response function is obtained as follows. We apply an
external field with the Gaussian profile described as

$\displaystyle\bm{F}(t)$
$\displaystyle=\frac{\bm{F}_{0}}{\sqrt{2\pi\sigma^{2}}}\mathrm{exp}\left(-\frac{(t-t_{0})^{2}}{2\sigma^{2}}\right).$
(15)

Since we took the velocity gauge by which
$\bm{E}(t)=-\frac{\partial\bm{A}(t)}{\partial t}$ holds, the Gaussian electric
field $\bm{E}(t)$ is incorporated by the vector potential $\bm{A}(t)$ given by

$\displaystyle\bm{A}(t)=\frac{\bm{E}_{0}}{2}\left[\mathrm{erf}\left(-\frac{t_{0}}{\sqrt{2}\sigma}\right)-\mathrm{erf}\left(\frac{t-t_{0}}{\sqrt{2}\sigma}\right)\right],$
(16)

with the error function $\mathrm{erf}(\cdot)$. In this scheme, we can
calculate the linear response function $\chi_{\mathcal{O}F}(\omega)$ of the
physical quantity $\mathcal{O}(t)$ to the external field $F(t)$ by the Fourier
transformation

$\displaystyle\chi_{\mathcal{O}F}(\omega)$
$\displaystyle=\frac{\Delta\mathcal{O}(\omega)}{F(\omega)},$
$\displaystyle=\frac{1}{F_{0}}e^{\frac{\sigma^{2}\omega^{2}}{2}}e^{i\omega
t_{0}}\int_{0}^{\infty}\mathcal{O}(t)e^{-i\omega t}dt,$ (17)

where $\Delta\mathcal{O}(\omega)$ is the Fourier component of
$\Delta\mathcal{O}(t)=\mathcal{O}(t)-\mathcal{O}(0)$ and $F(\omega)$ is the
Fourier component of the external field
$F(t)=\frac{F_{0}}{\sqrt{2\pi\sigma^{2}}}\mathrm{exp}\left(-\frac{(t-t_{0})^{2}}{2\sigma^{2}}\right)$.
In this calculation, we use the parameter of the external field
$F_{0}=1.0\times 10^{-5},~{}t_{0}=0.2,~{}\sigma=0.03$. Based on this scheme,
we calculate the linear response functions to the light field, which will be
shown in the next section.

We set the $k$-mesh of the Brillouin zone to $512\times 512$ for the real-time
simulation and the $k$-mesh of the Brillouin zone to $1600\times 1600$ for the
cubic interpolation of the SPDM in equilibrium. In this calculation, we
implement the parallel computation for $k$-mesh by the message passing
interface (MPI).

##  III Symmetry Analysis

Table 1: Eigenvalue of symmetry-adapted bases of localized spins under a symmetry operation in magnetic point group $\mathcal{G}$. The subscript $-$ of $\psi_{\alpha}$ represents the clockwise (CW) rotation, and the index $+$ represents the counterclockwise (CCW). In the columns of $3_{z}$ and $\theta 2_{x}$, $+1$ means the basis does not change its sign, while $-1$ means that the operation flips the sign. $\xi^{-}=\mathrm{exp}(2\pi i/3)$ denotes the eigenvalue of the CW mode and $\xi^{+}=\mathrm{exp}(-2\pi i/3)$ is that of the CCW mode. Each basis function $\psi_{\alpha}$ is explicitly written in the column ‘Basis function’.  $\psi_{\alpha}$ | $3_{z}$ | $\theta 2_{x}$ | Basis function
---|---|---|---
$\mathrm{A}_{1}^{1z}$ | $+1$ | $+1$ | $\mathrm{S}_{4}^{z}$
$\mathrm{A}_{1}^{3z}$ | $+1$ | $+1$ | $\frac{1}{\sqrt{3}}(\mathrm{S}_{1}^{z}+\mathrm{S}_{2}^{z}+\mathrm{S}_{3}^{z})$
$\mathrm{A}_{1}^{3xy}$ | $+1$ | $+1$ | $\frac{1}{\sqrt{3}}\left[(\frac{\sqrt{3}}{2}\mathrm{S}_{1}^{x}-\frac{1}{2}\mathrm{S}_{1}^{y})+(-\frac{\sqrt{3}}{2}\mathrm{S}_{2}^{x}-\frac{1}{2}\mathrm{S}_{2}^{y})+\mathrm{S}_{3}^{y}\right]$
$\mathrm{A}_{2}^{3xy}$ | $+1$ | $-1$ | $\frac{1}{\sqrt{3}}\left[(-\frac{1}{2}\mathrm{S}_{1}^{x}-\frac{\sqrt{3}}{2}\mathrm{S}_{1}^{y})+(-\frac{1}{2}\mathrm{S}_{2}^{x}+\frac{\sqrt{3}}{2}\mathrm{S}_{2}^{y})+\mathrm{S}_{3}^{x}\right]$
$\mathrm{E}_{\pm}^{1xy}$ | $\xi^{\pm}$ | $-1$ | $\frac{1}{\sqrt{2}}(\mathrm{S}_{4}^{x}\pm i\mathrm{S}_{4}^{y})$
$\mathrm{E}_{\pm}^{3z}$ | $\xi^{\pm}$ | $+1$ | $\frac{1}{\sqrt{2}}\left[\left(-\frac{1}{\sqrt{6}}\mathrm{S}_{1}^{z}-\frac{1}{\sqrt{6}}\mathrm{S}_{2}^{z}+\frac{2}{\sqrt{6}}\mathrm{S}_{3}^{z}\right)\pm i\left(-\frac{1}{\sqrt{2}}\mathrm{S}_{1}^{z}+\frac{1}{\sqrt{2}}\mathrm{S}_{2}^{z}\right)\right]$
$\mathrm{E}_{\pm}^{3xy}$ | $\xi^{\pm}$ | $+1$ | $\frac{1}{\sqrt{2}}\left[\frac{1}{\sqrt{3}}\left\\{(-\frac{\sqrt{3}}{2}\mathrm{S}_{1}^{x}-\frac{1}{2}\mathrm{S}_{1}^{y})+(\frac{\sqrt{3}}{2}\mathrm{S}_{2}^{x}-\frac{1}{2}\mathrm{S}_{2}^{y})+\mathrm{S}_{3}^{y}\right\\}\pm\frac{i}{\sqrt{3}}\left\\{(-\frac{1}{2}\mathrm{S}_{1}^{x}+\frac{\sqrt{3}}{2}\mathrm{S}_{1}^{y})+(-\frac{1}{2}\mathrm{S}_{2}^{x}-\frac{\sqrt{3}}{2}\mathrm{S}_{2}^{y})+\mathrm{S}_{3}^{x}\right\\}\right]$
$\mathrm{E}_{\pm}^{3xy^{\prime}}$ | $\xi^{\pm}$ | $-1$ | $\frac{1}{\sqrt{2}}\left[\frac{1}{\sqrt{3}}(\mathrm{S}_{1}^{x}+\mathrm{S}_{2}^{x}+\mathrm{S}_{3}^{x})\pm\frac{i}{\sqrt{3}}(\mathrm{S}_{1}^{y}+\mathrm{S}_{2}^{y}+\mathrm{S}_{3}^{y})\right]$

In this section, we present the symmetry analysis of the localized spin
system. First, we classify the collective modes of the localized spins.
Second, we analyze the collective spin dynamics that are linearly coupled to
the external field based on the symmetry analysis.

###  III.1 Symmetry-adapted basis of localized spins

In this subsection, we analyze the collective modes of the localized spins in
the framework of the magnetic representation and decompose them into three
modes: the azimuthal symmetric mode (AS), clockwise mode (CW), and
counterclockwise mode (CCW).

The system belongs to the magnetic point group $\mathcal{G}=32^{\prime}$ and
has the following symmetry

$\displaystyle\mathcal{G}=\\{1,3_{z},3_{z}^{-1},\theta 2_{x},\theta
2_{u},\theta 2_{v}\\}.$ (18)

Here, $1,\theta,3_{a},2_{a}$ are the identity operator, time-reversal
operator, threefold rotation around the $a$ axis, and twofold rotation around
the $a$ axis, respectively. The $u$ and $v$ axes are generated by the action
$3_{z},3_{z}^{-1}$ on the $x$ axis, respectively. The representation $\Gamma$
is obtained for the basis spanned by the degrees of freedom about localized
spin vectors $\mathrm{S}_{a}^{\mu}$. The dimension of $\Gamma$ is 12 due to
the sublattice $(a=1,2,3,4)$ and spins $(\mu=x,y,z)$. The representation is
decomposed into the irreducible representations as

$\displaystyle\Gamma=3\mathrm{A}_{1}\oplus\mathrm{A}_{2}\oplus 4\mathrm{E}.$

There are three symmetry-adapted bases for the irreducible representation
$\mathrm{A}_{1}$, one for $\mathrm{A}_{2}$, and eight for $\mathrm{E}$.

We further classify the symmetry-adapted modes by using the three-fold
rotations. We define the AS mode by the mode with its eigenvalue $1$ under the
symmetry operation $3_{z}$, the CW mode as that with
$\xi^{-}=\mathrm{exp}(2\pi i/3)$, and the CCW mode as that with
$\xi^{+}=\mathrm{exp}(-2\pi i/3)$. The symmetry-adapted basis function is
given by the linear combination of the localized spins $\mathrm{S}_{a}^{\mu}$
as

$\displaystyle\psi_{\alpha}=\sum_{a\mu}c^{\mu}_{a}(\alpha)\mathrm{S}_{a}^{\mu},$
(19)

where $c_{a}^{\mu}(\alpha)$ is the coefficient of the symmetry-adapted basis
$\psi_{\alpha}$.

The dynamics of the localized spin system can be described by the symmetry-
adapted basis listed in TABLE 1. We also display the sketch of the symmetry-
adapted basis functions $\left\\{\psi_{\alpha}\right\\}$ in FIG. 2. The AS
bases, $\mathrm{A}_{1}^{1z}$, $\mathrm{A}_{1}^{3z}$, $\mathrm{A}_{1}^{3xy}$,
and $\mathrm{A}_{2}^{3xy}$, are invariant under the out-of-plane rotation and
thus are totally symmetric in the azimuthal plane. The modes for clockwise
rotation include $\mathrm{E}_{-}^{1xy}$, $\mathrm{E}_{-}^{3z}$,
$\mathrm{E}_{-}^{3xy}$, and $\mathrm{E}_{-}^{3xy^{\prime}}$, while the CCW
mode comprises $\mathrm{E}_{+}^{1xy}$, $\mathrm{E}_{+}^{3z}$,
$\mathrm{E}_{+}^{3xy}$, and $\mathrm{E}_{+}^{3xy^{\prime}}$. Thus, CW and CCW
modes carry the nontrivial transformation property under the out-of-plane
rotation. Note that the CCW basis with $\xi_{+}=\mathrm{exp}(-2\pi i/3)$ is
the complex conjugate of the CCW counterpart with $\xi_{-}=\mathrm{exp}(2\pi
i/3)$. For example, $\mathrm{E}_{+}^{1xy}$ is the complex conjugate of
$\mathrm{E}_{-}^{1xy}$. These collective modes are similar to those of the
skyrmionic system [24, 25].

In the following sections, we calculate the component of the symmetry-adapted
basis of the localized spin dynamics $\tilde{\psi}_{\alpha}(t)$ as

$\displaystyle\tilde{\psi}_{\alpha}(t)=\sum_{a\mu}\left(c^{\mu}_{a}(\alpha)\right)^{*}\mathrm{S}_{a}^{\mu}(t),$
(20)

where we project the localized spins $\bm{\mathrm{S}}_{a}^{\mu}(t)$ into the
symmetry-adapted basis $\psi_{\alpha}$.

###  III.2 Selection rule for optical excitation of spin dynamics

Table 2: Eigenvalue of external field of localized spins under symmetry operator of magnetic point group $\mathcal{G}$. The symmetry-adapted bases listed in the column of the basis function are linearly coupled to the external field in the row. | $3_{z}$ | $\theta 2_{x}$ | Basis function
---|---|---|---
$E^{\pm}$ | $\xi^{\pm}$ | $-1$ | $\mathrm{E}_{\pm}^{1xy},\mathrm{E}_{\pm}^{3z},\mathrm{E}_{\pm}^{3xy},\mathrm{E}_{\pm}^{3xy^{\prime}}$
$B^{\pm}$ | $\xi^{\pm}$ | $+1$ | $\mathrm{E}_{\pm}^{1xy},\mathrm{E}_{\pm}^{3z},\mathrm{E}_{\pm}^{3xy},\mathrm{E}_{\pm}^{3xy^{\prime}}$
$B^{z}$ | $+1$ | $+1$ | $\mathrm{A}_{1}^{1z},\mathrm{A}_{1}^{3z},\mathrm{A}_{1}^{3xy},\mathrm{A}_{2}^{3xy}$
Figure 2:  Sketch of the symmetry-adapted basis of localized spins tabulated
in TABLE 2. The spin moments are colored as done in FIG. 1.

In this subsection, we analyze the symmetry-adapted bases of the localized
spins that are linearly coupled to the external field based on the symmetry.
The circularly polarized external field is defined as

$\displaystyle\bm{F}_{\pm}(t)=\frac{1}{\sqrt{2}}(1,\pm
i,0)\,\mathrm{exp}(i\omega t),$

where the minus sign ($-$) represents the CW mode and the plus ($+$)
represents the CCW mode. Similarly to TABLE 1, the external fields are
classified in terms of irreducible representations and eigenvalues for the
operation $3_{z}$. The result is shown in TABLE 2. For instance, the
eigenvalue of the clockwise rotatory field $\bm{F}_{-}(t)$ is $\xi^{-}$ for
the operation $3_{z}$, equal to that of the CW mode.

The symmetry indicates that we can excite the CW mode by the clockwise
rotatory field $\bm{F}_{-}(t)$ and the CCW mode by the counterclockwise
rotatory field $\bm{F}_{+}(t)$. The linearly polarized field $F_{x}$ contains
the clockwise and the counterclockwise rotation fields and can excite both the
CW and CCW modes. On the other hand, the eigenvalue of the magnetic field
$B_{z}$ is $1$ for $3_{z}$, which is the same as that of the AS mode. Thus, we
can excite the AS mode by the magnetic field $B_{z}$ along the $z$ axis. In
the following parts, we mainly delve into the linear responses to in-plane
electric fields, and the AS modes play minor roles.

##  IV Linear Response Functions

In this section, we present the results of the real-time simulation of the
spin-charge coupled system. First, we calculate the magnetic susceptibility of
the symmetry-adapted bases of the localized spin system in Sec. IV.1. Second,
we calculate the linear optical conductivity of the electron system in Sec.
IV.2. Third, we calculate the electromagnetic susceptibility of the localized
spin system in Sec. IV.3. We also quantitatively evaluate the electromagnetic
susceptibility of the collective modes.

###  IV.1 Magnetic Susceptibility

Figure 3: (a) Linear magnetic susceptibility of the AS mode to the out-of-
plane magnetic field. (b) Linear magnetic susceptibility of the CW mode to the
magnetic field along the $x$ axis. (c) Linear magnetic susceptibility of the
CCW mode to the in-plane magnetic field.

In this subsection, we show the magnetic susceptibility of the symmetry-
adapted bases of the localized spins to corroborate the collective spin
excitations in light of the frequency dependence. We calculate the magnetic
susceptibility $\chi_{\alpha\mu}(\omega)$ written in the symmetry-adapted
basis as

$\displaystyle\Delta\tilde{\psi}_{\alpha}(\omega)=\chi_{\alpha\mu}(\omega)B_{\mu}(\omega),$
(21)

where we defined the Fourier component of
$\Delta\tilde{\psi}_{\alpha}(t)=\tilde{\psi}_{\alpha}(t)-\tilde{\psi}_{\alpha}(0)$,
representing the modulation of localized spins in the symmetry-adapted
representation. The response function $\chi_{\alpha\mu}(\omega)$ is obtained
by taking the Gaussian magnetic field $\bm{F}(t)=\bm{B}(t)$ [Eq. (15)], which
is coupled to the localized spins. The perturbation is expressed by the
Hamiltonian

$H_{\text{ext}}(t)=-\sum_{a}\bm{B}(t)\cdot\bm{\mathrm{S}}_{a}.$ (22)

The perturbation takes place without breaking the translation symmetry, and
the time-evolution can be traced in parallel between different momentum as in
the case of the uniform vector potential [Eq. (6)].

We plot the magnetic susceptibility of the symmetry-adapted basis belonging to
the different symmetry in FIG. 3. In FIG. 3 (a), the magnetic susceptibility
$\chi_{\alpha z}$ with
$\alpha=\mathrm{A}_{1}^{1z},\mathrm{A}_{1}^{3z},\mathrm{A}_{1}^{3xy},\mathrm{A}_{2}^{3xy}$
(AS modes) are shown. Consistent with the symmetry analysis in Sec. III.2, the
out-of-plane magnetic field is coupled to the AS modes but makes no
contribution to the CW and CCW modes. On the other hand, the CW and CCW modes
respond to the in-plane magnetic fields as observed in the plots of
$\chi_{\alpha x}$ of FIG. 3 (b,c).

As a result, the spectrum of magnetic susceptibility is determined by three
types of collective modes and consists of four peaks, that is one for the AS
modes, two for the CW, and one for the CCW. The number of the collective modes
is consistent with the analysis of magnon spectrum [30], and the lowest-energy
($\omega\simeq 0.09$) peak comes from the magnon band whose gap energy at
$\bm{k}=\bm{0}$ is due to the magnetic anisotropy $K_{z}$ in Eq. (1). We note
that only the $\bm{k}=\bm{0}$ spin excitation is present since the system is
perturbed by the uniform light field parametrized by $\bm{E}(t)$.

###  IV.2 Optical Conductivity

In this subsection, we investigate the linear optical conductivity. The linear
optical conductivity is written by

$\displaystyle j_{\mu}(\omega)=\sigma_{\mu\nu}(\omega)E_{\nu}(\omega),$ (23)

where $\mu,\nu=x,y$ since the electronic Hamiltonian is built onto the two-
dimensional system. $j_{\mu}(\omega)$ is the Fourier component of the electric
current of Eq. (14).

Figure 4: (a) Linear optical conductivity of the system. The blue solid line
and the green dashed line indicate the longitudinal optical conductivity with
and without updating the spin configurations, respectively. The orange solid
line and the purple dashed line indicate the transverse optical conductivity
with and without updating the spin configurations, respectively. The black
dashed line indicates the band gap frequency. (b) Longitudinal optical
conductivity at log scale. The blue solid line and the orange dashed line
indicate the calculation with and without updating the spin configurations.

We plot the frequency dependence of the linear optical conductivity
$\sigma_{\mu\nu}(\omega)$ in FIG. 4. The optical gap is at $\omega\simeq 0.55$
[vertical dashed line in FIG. 4(a)], indicating the insulating state of the
adopted model. In FIG. 4(a), we show the longitudinal optical conductivity
$\sigma_{xx}$ and the transverse optical conductivity $\sigma_{yx}$. The solid
lines (‘With LLG’) represent the optical conductivity incorporating the effect
of the localized spin dynamics following the LLG equation, while the dashed
lines (‘Without LLG’) represent the optical conductivity without the effect.
In FIG. 4(b), we show the longitudinal optical conductivity $\sigma_{xx}$ with
and without the LLG simulation in log scale.

There are two consequences of spin dynamics on the optical conductivity.
First, $\operatorname{Re}\sigma(\omega)$ with the LLG simulation shows the
three resonance peaks in the in-gap regime, which are absent in the
independent particle approximation. Due to the increased degrees of freedom of
the local spins, more optical conductivity peaks are observed in the in-gap
regime than in the studies [4, 6]. The peaks at $\omega\simeq 0.09,0.25$
correspond to CW modes as inferred from the magnetic susceptibility [FIG. 3
(b)], and the peak at $\omega\simeq 0.16$ is for CCW modes [see FIG. 3(c)].
The peak height for $\omega\simeq 0.16$ is much smaller than that at
$\omega\simeq 0.09,0.25$. The prominent difference can be understood by the
electromagnetic susceptibility corroborated in Sec. IV.3.

Second, the spin dynamics leads to the modification of the optical
conductivity in the spectrum above the band gap as well as the in-gap optical
excitations. This modification is also carried by the collective spin
dynamics, while the localized spins are driven by the electronic spin
excitations responding to the irradiating light. It is noteworthy that the in-
gap and above-gap optical excitations are in the trade-off relationship due to
the sum rule
$\int\operatorname{Re}\sigma_{\mu\nu}(\omega)d\omega=\mathrm{const.}$ In the
spin-charge coupled system with the AIAO spin texture, the collective modes
belonging to the different symmetry affect the optical conductivity.
Therefore, we can selectively excite the collective modes which influence
electronic transport.

###  IV.3 Electromagnetic Susceptibility

Our focus is on the richer collective spin excitations originating from the
AIAO state and on their impact on the electric response, and thus let us make
more detailed investigations of the in-gap optical excitations by considering
the electromagnetic susceptibility. The electromagnetic susceptibility
$\kappa_{\alpha\mu}(\omega)$ of the symmetry-adapted form of localized spin
excitation $\psi_{\alpha}$ as

$\displaystyle\Delta\tilde{\psi}_{\alpha}(\omega)=\kappa_{\alpha\mu}(\omega)E_{\mu}(\omega).$
(24)

We plot the spectrum of the electromagnetic susceptibility $\kappa_{\alpha x}$
with the electric field $\bm{E}\parallel\hat{x}$ in FIG. 5. Owing to the
selection rule tabulated in TABLE III.2, only the CW and CCW modes contribute
to the electromagnetic susceptibility $\kappa_{\alpha x}$. Then, we show the
electromagnetic susceptibility relevant to the CW modes [FIG. 5(a)] and CCW
modes [FIG. 5(b)]. The electromagnetic susceptibility of the CCW mode is much
smaller than that of the CW mode.

The in-gap peaks of the optical conductivity are caused by resonant excitation
of collective spin dynamics, and thus peak amplitude implies the electrical
activity of each collective excitation. For instance, the electrical activity
of the CCW mode is small as observed in the electromagnetic susceptibility
plotted in FIG. 5(b), leading to the negligibly small peak in the optical
conductivity spectrum [FIG. 4(b)]. More specifically, the mechanism of
electric-active collective spins dynamics is as follows: first, the electric
field coupled to the electron system stimulates the magnetization response of
electrons’ spins through the spin-orbit coupling [third term in Eq. (1)].
Then, the perturbed spin moment of electrons induces the dynamics of the
localized spins as a consequence of the exchange interaction between them
[second term in Eq. (1)]. The resultant coupling between the electric field
and the localized spins indicates the electromagnetic response denoted by Eq.
(24).

Figure 5:  Linear electromagnetic susceptibility $\kappa_{\alpha x}$. The
symmetry-adapted basis $\psi_{\alpha}$ is taken from (a) the CW modes and from
(b) the CCW modes.

The mechanism can be formulated by decomposing the electromagnetic
susceptibility into two parts as

$\displaystyle\kappa_{\alpha x}(\omega)$
$\displaystyle=\sum_{\beta}\chi^{\bm{\mathrm{S}}}_{\alpha\beta}(\omega)\kappa^{\bm{\sigma}}_{\beta
x}(\omega),$ (25) $\displaystyle\equiv\sum_{\beta}\kappa_{\alpha\beta
x}(\omega).$ (26)

$\kappa_{\beta x}^{\bm{\sigma}}(\omega)$ is the bare electromagnetic
susceptibility indicating the correlation between the electric field as

$\sum_{\mu}c_{a}^{\mu}(\alpha)\Delta\Braket{\sigma_{a}^{\mu}}(\omega)=\kappa_{\alpha\nu}^{\bm{\sigma}}(\omega)E_{\nu}(\omega),$
(27)

where we define the Fourier component of
$\Delta\Braket{\sigma_{a}^{\mu}}(t)=\Braket{\sigma_{a}^{\mu}}(t)-\Braket{\sigma_{a}^{\mu}}(0)$,
representing the modulation of itinerant spins in the symmetry-adapted
representation. This susceptibility does not include the dynamical effects of
the localized spins, so we calculate $\kappa_{\alpha x}^{\bm{\sigma}}(\omega)$
by the real-time simulation without updating the local spin moments.

We also define $\chi_{\alpha\beta}^{\bm{\mathrm{S}}}(\omega)$, that is the
magnetic susceptibility denoting the correlation between the symmetry-adapted
modes ($\psi_{\alpha}$, $\psi_{\beta}$). Based on the symmetry,
$\chi_{\alpha\beta}^{\bm{\mathrm{S}}}(\omega)$ is finite when the symmetry-
adapted bases $\alpha,\beta$ belong to the same mode. The response function
$\chi_{\alpha\beta}^{\bm{\mathrm{S}}}(\omega)$ is obtained by taking the
Gaussian magnetic field, which is proportional to the symmetry-adapted basis
$\psi_{\beta}$. The perturbation is expressed by the Hamiltonian

$H_{\mathrm{ext}}^{\beta}(t)=-J\sum_{a\mu}\frac{B_{0}}{\sqrt{2\pi\sigma^{2}}}\mathrm{exp}\left(-\frac{(t-t_{0})^{2}}{2\sigma^{2}}\right)c_{a}^{\mu}(\beta)\mathrm{S}_{a}^{\mu},$
(28)

where $B_{0}$ is the amplitude of the Gaussian magnetic field.

Combining $\chi_{\alpha\beta}^{\bm{\mathrm{S}}}(\omega)$ and $\kappa_{\beta
x}^{\bm{\sigma}}(\omega)$, we can calculate the decomposed component
$\kappa_{\alpha\beta x}(\omega)$. As a result, the electromagnetic
susceptibility $\kappa_{\alpha\mu}(\omega)$ of the collective mode
$\psi_{\alpha}$ is determined by the interference of
$\kappa_{\alpha\beta\mu}(\omega)$ with respect to the mediating modes labeled
by $\beta$.

Figure 6: (a) The decomposition of the electromagnetic susceptibility of the
collective mode $\psi_{\alpha}=\mathrm{E}_{-}^{1xy}$. (b) The decomposition of
the electromagnetic susceptibility of the collective mode
$\psi_{\alpha}=\mathrm{E}_{+}^{3xy^{\prime}}$.

Taking the CW $\mathrm{E}_{-}^{1xy}$ and CCW modes
$\mathrm{E}_{+}^{3xy^{\prime}}$ for $\psi_{\alpha}$ in Eq. (26), we plot the
decomposed electromagnetic susceptibility in FIG. 6. In the case of
$\psi_{\alpha}=\mathrm{E}_{-}^{1xy}$ for example, the localized-spin response
to the electric field is mediated by the fluctuations denoted by
$\sum_{\mu}c_{a}^{\mu}(\beta)\Braket{\sigma_{a}^{\mu}}$ which are symmetry-
adapted to the CW modes ($\mathrm{E}_{-}^{1xy}$, $\mathrm{E}_{-}^{3z}$,
$\mathrm{E}_{-}^{3xy}$, and $\mathrm{E}_{-}^{3xy^{\prime}}$). Similarly, the
electromagnetic susceptibility of the collective mode
$\mathrm{E}_{+}^{3xy^{\prime}}$ is mediated by those sharing the same symmetry
with the CCW modes of localized spins.

As for the electromagnetic susceptibility of the CW mode, the interference of
the mediating spin fluctuations results in a sizable electromagnetic
susceptibility for $\psi_{\alpha}=\mathrm{E}_{-}^{1xy}$. On the other hand,
while each $\kappa_{\alpha\beta x}$ makes a significant contribution,
destructive interference is observed in the case of the CCW mode, giving rise
to a much smaller electromagnetic susceptibility in total. These results imply
that the CCW modes are less electrically active when compared to the CW modes
and hence offer negligible optical conductivity [peak around $\omega\simeq
0.16$ in FIG. 4(b)]. The constructive and destructive interference effects are
features unique to the system hosting the complex spin texture, and their
identification has been accomplished by the developed efficient computational
method for the real-time evolution of the spin-charge coupled system.

##  V Summary

In this study, we developed an efficient real-time simulation of the spin-
charge coupled system and applied this method to the two-dimensional system
with a complex spin texture. We investigated how the collective excitation of
the localized spins influences the optical response of the electronic system.
First, we identified the symmetry-adapted bases of the localized spins and
further classified them with the circular polarization. This allowed us to
clarify the activity of the collective modes of the localized spins to
external stimuli. Next, we calculated the linear response functions of the
spin-charge coupled system. The magnetic susceptibility is conveniently
decomposed by the symmetry-adapted bases and is determined by the resonant
dynamics of the localized spins when the frequency of magnetic fields is in
the in-gap regime. The observed collective modes have features similar to
those of the skyrmionic system, implying the AIAO state is the minimum unit of
the magnetic skyrmion.

The impact of the real-time simulation has been demonstrated by the
calculations of the linear optical conductivity including the effect of the
localized-spins dynamics. The electrical activity of the collective spin
motion results in the modification of the spectrum of optical conductivity
such as in-gap peak structures. The collective modes differ in their influence
on the optical conductivity; the rotatory modes of the localized spins affect
the optical conductivity, and the AS mode of the localized spins does not
affect the optical conductivity in the two-dimensional system. Finally, we
discussed electromagnetic susceptibility, paying special attention to spin
dynamics induced by the light whose frequency is below the band gap. The
electromagnetic susceptibilities for the CW and CCW modes are elucidated by
the decomposition into the bare electromagnetic susceptibility and magnetic
susceptibility. As a result, the interference effect plays a key role in
determining the contrasting electrical activity of the CW and CCW modes.

To conclude, we systematically investigated the spin-charge coupling by real-
time simulations of the system with the AIAO magnetic order, which is
considered the minimal unit of magnetic skyrmion. The diverse collective
excitations arising from the magnetic structure lead to multiple peak
structures in the optical conductivity spectrum and exhibit different
electrical activities due to interference effects between spin fluctuations.
The essential ingredient is the optically active collective spin dynamics in
the gapped system. Therefore, searching for a real material should be in the
realm of inversion broken semiconductors or insulators with optically active
magnons, such as $\mathrm{GaV_{4}S_{8}}$ [34, 35] and $\mathrm{Cu_{2}OSeO}$
[36, 37].

## acknowledgement

This work is supported by Grant-in-Aid for Scientific Research from JSPS,
KAKENHI Grant No. JP23K13058 (H.W.), No. 22H00290(T.N.), No. 24K00581(T.N.),
No. 21H04990 (R.A.), No. JP21H04437 (R.A.), No. 19H05825 (R.A.), JST-PRESTO
No. JPMJPR20L7 (T.N.), JST-CREST No. JPMJCR18T3, No. JPMJCR23O4(R.A.), JST-
ASPIRE No. JPMJAP2317 (R.A.), JST-Mirai No. JPMJMI20A1 (R.A.). K.H. was
supported by the Program for Leading Graduate Schools (MERIT-WINGS).

## References

  * Sotome  _et al._ [2021] M. Sotome, M. Nakamura, T. Morimoto, Y. Zhang, G.-Y. Guo, M. Kawasaki, N. Nagaosa, Y. Tokura, and N. Ogawa, Physical Review B 103, L241111 (2021).
  * Morimoto and Nagaosa [2016] T. Morimoto and N. Nagaosa, Physical Review B 94, 035117 (2016).
  * Okamura  _et al._ [2022] Y. Okamura, T. Morimoto, N. Ogawa, Y. Kaneko, G.-Y. Guo, M. Nakamura, M. Kawasaki, N. Nagaosa, Y. Tokura, and Y. Takahashi, Proceedings of the National Academy of Sciences 119, e2122313119 (2022).
  * Iguchi  _et al._ [2024] J. Iguchi, H. Watanabe, Y. Murakami, T. Nomoto, and R. Arita, Physical Review B 109, 064407 (2024).
  * Murakami  _et al._ [2020] Y. Murakami, D. Golež, T. Kaneko, A. Koga, A. J. Millis, and P. Werner, Physical Review B 101, 195118 (2020).
  * Kaneko  _et al._ [2021] T. Kaneko, Z. Sun, Y. Murakami, D. Golež, and A. J. Millis, Physical Review Letters 127, 127402 (2021).
  * Ono and Ishihara [2021] A. Ono and S. Ishihara, npj Computational Materials 7, 171 (2021).
  * Ohgushi  _et al._ [2000] K. Ohgushi, S. Murakami, and N. Nagaosa, Physical Review B 62, R6065 (2000).
  * Martin and Batista [2008] I. Martin and C. D. Batista, Physical Review Letters 101, 156402 (2008).
  * Nagaosa and Tokura [2012] N. Nagaosa and Y. Tokura, Physica Scripta 2012, 014020 (2012).
  * Nagaosa and Tokura [2013] N. Nagaosa and Y. Tokura, Nature Nanotechnology 8, 899 (2013).
  * Mühlbauer  _et al._ [2009] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Science 323, 915 (2009).
  * Kurumaji  _et al._ [2019] T. Kurumaji, T. Nakajima, M. Hirschberger, A. Kikkawa, Y. Yamasaki, H. Sagayama, H. Nakao, Y. Taguchi, T. hisa Arima, and Y. Tokura, Science 365, 914 (2019).
  * Hirschberger  _et al._ [2019] M. Hirschberger, T. Nakajima, S. Gao, L. Peng, A. Kikkawa, T. Kurumaji, M. Kriener, Y. Yamasaki, H. Sagayama, H. Nakao, K. Ohishi, K. Kakurai, Y. Taguchi, X. Yu, T.-h. Arima, and Y. Tokura, Nature Communications 10, 5831 (2019).
  * Akagi and Motome [2010] Y. Akagi and Y. Motome, Journal of the Physical Society of Japan 79, 083711 (2010).
  * Hamamoto  _et al._ [2015] K. Hamamoto, M. Ezawa, and N. Nagaosa, Physical Review B 92, 115417 (2015).
  * Matsui  _et al._ [2021] A. Matsui, T. Nomoto, and R. Arita, Physical Review B 104, 174432 (2021).
  * Feng  _et al._ [2020] W. Feng, J.-P. Hanke, X. Zhou, G.-Y. Guo, S. Blügel, Y. Mokrousov, and Y. Yao, Nature Communications 11, 118 (2020).
  * Zhou  _et al._ [2016] J. Zhou, Q.-F. Liang, H. Weng, Y. B. Chen, S.-H. Yao, Y.-F. Chen, J. Dong, and G.-Y. Guo, Physical Review Letters 116, 256601 (2016).
  * Park  _et al._ [2023] P. Park, W. Cho, C. Kim, Y. An, Y.-G. Kang, M. Avdeev, R. Sibille, K. Iida, R. Kajimoto, K. H. Lee, W. Ju, E.-J. Cho, H.-J. Noh, M. J. Han, S.-S. Zhang, C. D. Batista, and J.-G. Park, Nature Communications 14, 8346 (2023).
  * Takagi  _et al._ [2023] H. Takagi, R. Takagi, S. Minami, T. Nomoto, K. Ohishi, M.-T. Suzuki, Y. Yanagi, M. Hirayama, N. D. Khanh, K. Karube, H. Saito, D. Hashizume, R. Kiyanagi, Y. Tokura, R. Arita, T. Nakajima, and S. Seki, Nature Physics 19, 961 (2023).
  * Ghimire  _et al._ [2018] N. J. Ghimire, A. S. Botana, J. S. Jiang, J. Zhang, Y.-S. Chen, and J. F. Mitchell, Nature Communications 9, 3280 (2018).
  * Khanh  _et al._ [2020] N. D. Khanh, T. Nakajima, X. Yu, S. Gao, K. Shibata, M. Hirschberger, Y. Yamasaki, H. Sagayama, H. Nakao, L. Peng, K. Nakajima, R. Takagi, T.-h. Arima, Y. Tokura, and S. Seki, Nature Nanotechnology 15, 444 (2020).
  * Mochizuki [2012] M. Mochizuki, Physical Review Letters 108, 017601 (2012).
  * Mruczkiewicz  _et al._ [2017] M. Mruczkiewicz, M. Krawczyk, and K. Y. Guslienko, Physical Review B 95, 094414 (2017).
  * Mruczkiewicz  _et al._ [2016] M. Mruczkiewicz, P. Gruszecki, M. Zelent, and M. Krawczyk, Physical Review B 93, 174429 (2016).
  * Seki  _et al._ [2020] S. Seki, M. Garst, J. Waizner, R. Takagi, N. D. Khanh, Y. Okamura, K. Kondou, F. Kagawa, Y. Otani, and Y. Tokura, Nature Communications 11, 256 (2020).
  * Onose  _et al._ [2012] Y. Onose, Y. Okamura, S. Seki, S. Ishiwata, and Y. Tokura, Physical Review Letters 109, 037603 (2012).
  * Okamura  _et al._ [2013] Y. Okamura, F. Kagawa, M. Mochizuki, M. Kubota, S. Seki, S. Ishiwata, M. Kawasaki, Y. Onose, and Y. Tokura, Nature Communications 4, 2391 (2013).
  * Akagi  _et al._ [2013] Y. Akagi, M. Udagawa, and Y. Motome, Journal of the Physical Society of Japan 82, 123709 (2013).
  * Yue and Gaarde [2022] L. Yue and M. B. Gaarde, Journal of the Optical Society of America B 39, 535 (2022).
  * Ventura  _et al._ [2017] G. B. Ventura, D. J. Passos, J. M. B. L. dos Santos, J. M. V. P. Lopes, and N. M. R. Peres, Physical Review B 96, 035431 (2017).
  * Murakami and Schüler [2022] Y. Murakami and M. Schüler, Physical Review B 106, 035204 (2022).
  * Padmanabhan  _et al._ [2019] P. Padmanabhan, F. Sekiguchi, R. B. Versteeg, E. Slivina, V. Tsurkan, S. Bordács, I. Kézsmárki, and P. H. M. van Loosdrecht, Physical Review Letters 122, 107203 (2019).
  * Kézsmárki  _et al._ [2015] I. Kézsmárki, S. Bordács, P. Milde, E. Neuber, L. M. Eng, J. S. White, H. M. Rønnow, C. D. Dewhurst, M. Mochizuki, K. Yanai, H. Nakamura, D. Ehlers, V. Tsurkan, and A. Loidl, Nature Materials 14, 1116 (2015).
  * Seki  _et al._ [2012] S. Seki, X. Z. Yu, S. Ishiwata, and Y. Tokura, Science 336, 198 (2012).
  * White  _et al._ [2014] J. S. White, K. Prša, P. Huang, A. A. Omrani, I. Živković, M. Bartkowiak, H. Berger, A. Magrez, J. L. Gavilano, G. Nagy, J. Zang, and H. M. Rønnow, Phys. Rev. Lett. 113, 107203 (2014).


11institutetext: Instytut Astronomiczny, Uniwersytet Wrocławski, ul. Kopernika
11, 51-622 Wrocław, Poland
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Effect of NLTE model atmospheres on photometric amplitudes
and phases of early B-type pulsating stars

J. Daszyńska-Daszkiewicz  W. Szewczuk

(Received …; accepted …)

###### Abstract

Context. Amplitudes and phases of the light variation of a pulsating star in
various photometric passbands contain information about geometry of observed
modes. Because oscillation spectra of early B-type main sequence stars do not
exhibit regular patterns, these observables are very often the only ones from
which mode identification can be derived. Moreover, these data can yield
valuable constraints on mean stellar parameters, subphotospheric convection,
microphysics and atmospheres.

Aims. We study all possible sources of inaccuracy in theoretical values of the
photometric observables, i.e. amplitude ratios and phase differences, of early
B-type main sequence pulsators. Here, we discuss effects of parameters coming
from both models of stellar atmospheres and linear nonadiabatic theory of
stellar pulsation. In particular, we evaluate for the first time the effect of
the departure from the LTE approximation.

Methods. The photometric amplitudes and phases are calculated from the semi-
analytical formula for the light variation. The atmospheric input comes from
line-blanketed, LTE and NLTE plane-parallel, hydrostatic models. To compute
the limb-darkening coefficients for NLTE models, we use the Least-Square
Method taking into account the accuracy of the flux conservation. The linear
nonadiabatic stellar pulsations were computed by means of Dziembowski code. We
consider the OPAL and OP opacity tables and two determinations of the solar
chemical mixture: GN93 and AGSS09.

Results. We present effects of NLTE atmospheres, chemical composition and
opacities on theoretical values of the photometric observables of early B-type
pulsators. To this end, we compute tables with the passband fluxes, flux
derivatives over effective temperature and gravity as well as the non-linear
limb-darkening coefficients in 12 most often used passbands, i.e. in the
Strömgern system, $uvby$, and in the Johnson-Cousins-Glass system, $UBVRIJHK$.
We make these tables public available at the Wrocław HELAS Web page ††thanks:
http://helas.astro.uni.wroc.pl.

Conclusions. In the case of radial modes, effects of uncertainties in model
atmospheres are much smaller than those coming from the pulsation theory. In
turn, for nonradial modes effects of NLTE become more important and they are
most significant for dipole modes. Therefore, if one wants to construct an
accurate seismic model and gain more constraints, all inaccuracies in the
photometric observables have to be allowed for.

###### Key Words.:

Stars: early-type – Stars: oscillations – Stars: atmospheres

##  1 Introduction

Data on the photometric and spectroscopic variations of a pulsating star bring
information on frequencies of excited modes and their geometry. The latter
property is particularly important if oscillation spectra are sparse and lack
equidistant patterns. This is the case for main sequence pulsators of early B
spectral type, i.e. $\beta$ Cephei variables. Additionally, using data on the
photometric amplitudes and phases one can constrain mean stellar parameters,
test efficiency of convection transport in the case of cooler pulsators, like
$\delta$ Scuti stars, (Daszyńska-Daszkiewicz, Dziembowski, Pamyatnykh 2003) or
test opacity data in the case of B-type pulsators (Daszyńska-Daszkiewicz,
Dziembowski, Pamyatnykh 2005).

Ones of the most popular tools to identify a pulsation mode are the amplitude
ratios and phase differences in various photometric passbands. In the zero-
rotation approximation, these observables depend on the mode degree, $\ell$,
but are independent of the azimuthal order, $m$, the intrinsic mode amplitude,
$\varepsilon$, and the inclination angle, $i$.

The semi-analytical expression for the bolometric light variation was
formulated by Dziembowski (dziemb (1977)). Then, Balona & Stobie (balst
(1979)) and Stamford & Watson (stwa (1981)) expanded it for the light
variation in the photometric passbands. They showed that modes with different
values of $\ell$ are located in separated parts on the amplitude ratio $vs.$
phase difference diagrams. Subsequently, this method has been applied to
various types of pulsating stars by Watson (watson (1988)). Cugier,
Dziembowski & Pamyatnykh (cdp (1994)) improved the method by including
nonadiabatic effects in calculations for the $\beta$ Cephei stars. Effects of
rotation on photometric observables were studied by Daszyńska-Daszkiewicz et
al. (dd2002 (2002)) for close frequency modes and by Townsend (townsend2003
(2003)) and Daszyńska-Daszkiewicz, Dziembowski & Pamyatnykh (dd2007 (2007))
for long-period g-modes. To compute values of the photometric amplitudes and
phases we need two kinds of input. The first one comes from stellar model
atmospheres and the second one from computations of stellar pulsation. Both
data contain various sources of uncertainties.

The goal of this paper is to examine for the first time an influence of NLTE
effects on theoretical values of the photometric amplitude ratios and phase
differences for early B-type main sequence pulsators. We check also effects of
metallicity and microturbulent velocity in the atmosphere. Subsequently, we
compare these effects with uncertainties coming from the linear nonadiabatic
theory of stellar pulsation. As an example, we consider main sequence models
with a mass of 10 $M_{\odot}$ and low degree modes with $\ell$=0, 1, 2. Here,
we neglect all effects of rotation on pulsation.

The structure of the paper is the following. In Section 2, we recall the
linear formula for the pulsation complex amplitude in a photometric band.
Section 3 contains description of NLTE model atmospheres and results of our
computations of the band fluxes, corresponding flux derivatives over effective
temperature and gravity, and the nonlinear limb darkening coefficients, in 12
photometric passbands: $uvbyUBVRIJHK$. Tables with these data can be
downloaded from the Wrocław HELAS web page. Moreover, we study effects of
temperature, gravity, NLTE, atmospheric metallicity and microturbulent
velocity on the above mentioned quantities. How these atmospheric
uncertainties translate into the pulsation photometric observables, i.e.
amplitude ratios and phase differences, is discussed in Section 4.
Inaccuracies in the photometric observables connected with the linear
nonadiabatic theory of stellar pulsation are presented in Section 5. We end
with conclusions in Section 6.

##  2 Light variation due to stellar pulsation

Stellar pulsations cause changes of temperature, normal to the surface element
and pressure. If pulsations are linear and all effects of rotation on
pulsation can be ignored, then the total amplitude of the light variation in
the passband $x$ can be written in the following complex form (Daszyńska-
Daszkiewicz et al. dd2002 (2002)):

$\mathcal{A}_{x}(i)=-1.086\varepsilon
Y_{\ell}^{m}(i,0)b_{\ell}^{x}(D_{1,\ell}^{x}f+D_{2,\ell}+D_{3,\ell}^{x}),$ (1)

where $\varepsilon$ is the intrinsic mode amplitude, $Y_{\ell}^{m}$ – the
spherical harmonic and $i$ – the inclination angle. The amplitude is given by
$abs(\mathcal{A}_{x})$ and phase by $arg(\mathcal{A}_{x})$. The
$D_{1,\ell}^{x}\cdot f$ product stands for temperature changes, where

$D_{1,\ell}^{x}=\frac{1}{4}\frac{\partial\log(\mathcal{F}_{x}|b_{\ell}^{x}|)}{\partial\log
T_{\rm eff}}.$ (2)

${\cal F}_{x}$ is the flux in the passband $x$ and $f$ is the nonadiabatic
complex parameter describing the amplitude of the radiative flux perturbation
to the radial displacement at the photosphere level

$\frac{\delta{\cal F}_{\rm bol}}{{\cal F}_{\rm bol}}={\rm Re}\\{\varepsilon
fY_{\ell}^{m}(\theta,\varphi){\rm e}^{-{\rm i}\omega t}\\}.$ (3)

Geometrical term, $D_{2,\ell}$, is given by

$D_{2,\ell}=(2+\ell)(1-\ell),$ (4)

and the pressure term, $D_{3,\ell}^{x}$, by

$D_{3,\ell}^{x}=-\left(2+\frac{\omega^{2}R^{3}}{GM}\right)\frac{\partial\log(\mathcal{F}_{x}|b_{\ell}^{x}|)}{\partial\log
g}.$ (5)

$b_{\ell}^{x}$ is the disc averaging factor defined by

$b_{\ell}^{x}=\int_{0}^{1}h_{x}(\mu)\mu P_{\ell}(\mu)d\mu$ (6)

where $h_{x}(\mu)$ is the limb darkening law, $P_{\ell}$ is the Legendre
polynomial and $\mu$ is a cosine of the angle between a line of sight and the
emergent intensity. Remaining parameters have their usual meaning.

In the above expressions, we can distinguish two sorts of input parameters
needed to compute theoretical values of the photometric amplitudes and phases.
The first input is derived from models of stellar atmospheres and these are
the flux derivatives over effective temperature and gravity (Eq. 2 and 5), as
well as limb-darkening and its derivatives (Eq. 2, 5 and 6). The second input
comes from the nonadiabatic theory of stellar pulsation and this is the
$f$-parameter (Eq. 1 and 3).

In this paper, we used the Warsaw-New Jersey evolutionary code and the linear
nonadiabatic pulsation code of Dziembowski (1977). We considered opacity
tables from OPAL (Iglesias & Rogers 1996) and OP (Seaton 2005) projects, and
two determinations of the solar chemical composition: GN93 by Grevesse & Noels
1993 and AGSS09, a recent one by Asplund et al. 2009. As for models of stellar
atmospheres, we considered Kurucz models (Kurucz kurucz (2004)) and TLUSTY
models (Lanz & Hubeny lanzhubeny (2007)).

##  3 NLTE line-blanketed model atmospheres

The most widely used models of stellar atmospheres are line blanketed, plane-
parallel, hydrostatic models of Kurucz (2004) computed within an approximation
of local thermodynamic equilibrium (LTE). However, in the case of atmospheres
of early B-type stars, effects of the departure from LTE and a proper
treatment of line opacity become important. Escalation of the quality of
spectro-photometric observations calls for a need of high resolution and
accurate models of stellar atmospheres.

A grid of non-LTE (NLTE) model atmospheres were computed by Lanz & Hubeny
(lanzhubenyO (2003)) and Lanz & Hubeny (lanzhubeny (2007)). These are metal
line-blanketed, plane-parallel, hydrostatic model atmospheres of O-type stars
(OSTAR2002) and of early B-type stars (BSTAR2006), respectively. In their
computations, they adopted a solar chemical mixture by Grevesse & Sauval
(grevese (1998)), helium to hydrogen abundance of He/H=0.1 by number and two
values of the microturbulent velocity, $\xi_{t}$=2 and 10 km/s. Lanz & Hubeny
(lanzhubenyO (2003), lanzhubeny (2007)) showed that in the case of OB stars, a
neglect of NLTE effects causes differences not only in spectral lines but also
in the continuum level. In the near ultraviolet (the Balmer continuum), the
LTE fluxes are up to $10\%$ higher than the NLTE counterparts. In turn, in the
far and extreme ultraviolet (the Lyman continuum), the LTE fluxes are lower
than the NLTE ones. For more details see Lanz & Hubeny (lanzhubenyO (2003),
lanzhubeny (2007)).

The grid of the OSTAR2002 models were computed for 12 values of effective
temperatures appropriate to O-type stars, i.e. between 27500 and 55000 K with
a step of 2500 K, 8 values of surface gravities from $\log g$=3.0 to 4.75 with
a step of 0.25 dex and one microturbulent velocity, $\xi_{t}$=10 km/s.
Moreover, 10 values of the atmospheric metallicity were considered,
$(Z/Z_{\odot})_{\rm atm}$=2, 1, 0.5, 0.2, 0.1, 1/30, 1/50, 1/100, 1/1000, 0.0,
where $Z$ is the metal abundance by mass and $Z_{\odot}$ is the solar value.

The grid of the BSTAR2006 models contains 16 values of effective temperatures
between 15000 and 30000 K and a step of 1000 K, 13 values of surface gravities
in the range from $\log g=1.75$ to 4.75 dex and a 0.25 dex step. In this case,
6 values of metallicity were considered, $(Z/Z_{\odot})_{\rm atm}$=2, 1, 0.5,
0.2, 0.1, 0.0, and two values of the microturbulent velocities, $\xi_{t}$=2
km/s and $\xi_{t}$=10 km/s. Models with $\xi_{t}$=10 km/s were computed for
two sets of chemical mixture but only for B-type supergiants ($\log g\leq
3.0$). The first mixture was the same as for $\xi_{t}$=2 km/s, i.e. the solar
composition, and the second one was enriched in helium and nitrogen, and
depleted in carbon.

The lower limit for gravity at a given effective temperature was determined
approximately by the Eddington limit. The TLUSTY code used by Hubeny and Lanz
(1995) becomes unstable near this limit.

###  3.1 The passband fluxes and their derivatives

We computed the fluxes, ${\cal F}_{x}$, for the BSTAR2006 NLTE models in 12
commonly used photometric passbands, i.e. in the Strömgren ($uvby$) and
Johnson-Cousins-Glass ($UBVRIJHK$) systems, according to the formula

${\cal F}_{x}=\frac{\int\limits^{\lambda_{1}}_{\lambda_{2}}{\cal
F}(\lambda)S(\lambda)\,d\lambda}{\int\limits^{\lambda_{1}}_{\lambda_{2}}S(\lambda)\,d\lambda},$
(7)

where $S(\lambda)$ is the response function of the passband $x$, adopted from
the Asiago Database on Photometric Systems (Moro & Munari moro (2000)). The
integral in Eq. (7) is computed in the wavelength range from $\lambda_{1}$ to
$\lambda_{2}$ where $S(\lambda)$ has non-zero values.

We considered the whole grid of effective temperature, $T_{\rm eff}$, gravity,
$\log g$, metallicity, $(Z/Z_{\odot})_{\rm atm}$, and microturbulent velocity,
$\xi_{t}$, for the BSTAR2006 models. Table 1 summarizes our results. The file
names are coded by values of $(Z/Z_{\odot})_{\rm atm}$ and $\xi_{t}$. We have
adopted the same designations as Lanz &Hubeny (lanzhubeny (2007)), i.e. BC,
BG, BL, BS, BT, BZ denote $(Z/Z_{\odot})_{\rm atm}$=2, 1, 0.5, 0.2, 0.1 and 0,
respectively, whereas v2 and v10 correspond to $\xi_{t}=2$ km/s and 10 km/s,
respectively. The index CN marks a model enriched in helium and nitrogen, and
depleted in carbon. Each file contains the following columns: line number,
effective temperature, $T_{\rm eff}$, logarithm of the surface gravity, $\log
g$, metallicity, $(Z/Z_{\odot})_{\rm atm}$, microturbulent velocity,
$\xi_{t}$, and the logarithmic flux, $\log{\cal F}_{x}$, in the $uvbyUBVRIJHK$
passbands.

Table 1: Tables of the passband fluxes, their derivatives over effective
temperature and gravity, and non-linear limb darkening coefficients, for the
BSTAR2006 NLTE models in the $uvbyUBVRIJHK$ passbands.

fluxes | derivatives | limb darkening | range of $T_{\rm eff}$ [K] | range of $\log g$ [dex] | $(Z/Z_{\odot})_{\rm atm}$ |  $\xi_{t}$ [km/s] | mixture
---|---|---|---|---|---|---|---
|  | coefficients |  |  |  |  |
flux_BCv2 | der_BCv2 | LDC_BCv2 | 15000-30000 | 1.75-4.75 | 2.0 | 2 | GS98
flux_BCv10 | der_BCv10 | LDC_BCv10 | 15000-30000 | 1.75-3.00 | 2.0 | 10 | GS98
flux_BCv10CN | der_BCv10CN | LDC_BCv10CN | 15000-30000 | 1.75-3.00 | 2.0 | 10 | CN
flux_BGv2 | der_BGv2 | LDC_BGv2 | 15000-30000 | 1.75-4.75 | 1.0 | 2 | GS98
flux_BGv10 | der_BGv10 | LDC_BGv10 | 15000-30000 | 1.75-3.00 | 1.0 | 10 | GS98
flux_BGv10CN | der_BGv10CN | LDC_BGv10CN | 15000-30000 | 1.75-3.00 | 1.0 | 10 | CN
flux_BLv2 | der_BLv2 | LDC_BLv2 | 15000-30000 | 1.75-4.75 | 0.5 | 2 | GS98
flux_BLv10 | der_BLv10 | LDC_BLv10 | 15000-30000 | 1.75-3.00 | 0.5 | 10 | GS98
flux_BLv10CN | der_BLv10CN | LDC_BLv10CN | 15000-30000 | 1.75-3.00 | 0.5 | 10 | CN
flux_BSv2 | der_BSv2 | LDC_BSv2 | 15000-30000 | 1.75-4.75 | 0.2 | 2 | GS98
flux_BSv10 | der_BSv10 | LDC_BSv10 | 15000-30000 | 1.75-3.00 | 0.2 | 10 | GS98
flux_BSv10CN | der_BSv10CN | LDC_BSv10CN | 15000-30000 | 1.75-3.00 | 0.2 | 10 | CN
flux_BTv2 | der_BTv2 | LDC_BTv2 | 15000-30000 | 1.75-4.75 | 0.1 | 2 | GS98
flux_BTv10 | der_BTv10 | LDC_BTv10 | 15000-30000 | 1.75-3.00 | 0.1 | 10 | GS98
flux_BTv10CN | der_BTv10CN | LDC_BTv10CN | 15000-30000 | 1.75-3.00 | 0.1 | 10 | CN
flux_BZv2 | der_BZv2 | LDC_BZv2 | 15000-30000 | 1.75-4.75 | 0.0 | 2 | GS98
flux_BZv10 | der_BZv10 | LDC_BZv10 | 15000-30000 | 1.75-3.00 | 0.0 | 10 | GS98

111 GS98 - chemical mixture by Grevesse & Sauval (grevese (1998)), CN -
chemical mixture enriched in helium and nitrogen, and depleted in carbon.

The most important atmospheric parameters in the expression for the brightness
variation of a pulsating star are the flux derivatives over effective
temperature, $T_{\rm eff}$, and gravity, $\log g$. In a given photometric
pasband $x$, they are defined as

$\alpha_{T}^{x}=\frac{\partial\log{\cal F}_{x}}{\partial\log T_{\rm
eff}}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ {\rm
and}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\alpha_{g}^{x}=\frac{\partial\log{\cal F}_{x}}{\partial\log g}.$

Values of these flux derivatives depend not only on effective temperature and
gravity but also on the metallicity, microturbulent velocity and the departure
from LTE in the star’s atmosphere.

Similarly to the passband fluxes, the derivatives were computed in the whole
range of parameters of the BSTAR2006 models. Names of the derivative tables
are given in Table 1. Each file contains the following columns: line number,
effective temperature, $T_{\rm eff}$, logarithm of the surface gravity, $\log
g$, metallicity, $(Z/Z_{\odot})_{\rm atm}$, microturbulent velocity,
$\xi_{t}$, and the flux derivatives over effective temperature, $\alpha_{T}$,
and gravity, $\alpha_{g}$, in the $uvbyUBVRIJHK$ passbands.

In Fig. 1, we show the NLTE flux derivatives as a function of temperature for
two Strömgren passbands $uy$ and three values of gravity, $\log
g=3.5,4.0,4.5$. In the left panel, we plot the temperature derivative,
$\alpha_{T}$, and in the right one the gravity derivative, $\alpha_{g}$. We
assumed the solar metallicity, $(Z/Z_{\odot})_{\rm atm}$=1 and the
microturbulent velocity of $\xi_{t}$=2 km/s. As we can see the lower the
gravity the larger values of $\alpha_{T}$ and the lower values of
$\alpha_{g}$.

The wavelength dependence of the flux derivatives is presented in Fig. 2. We
considered the central wavelengths, $\lambda_{c}$, of the $UBVRIJHK$ passbands
and a model with $T_{\rm eff}=20000$ K and $\log g$ =4.0. The most steep
derivatives are in the ultraviolet filter and then the absolute values of
$\alpha_{T}$ and $\alpha_{g}$ decrease quite rapidly with $\lambda_{c}$. This
is because the early B-type stars emit the most amount of their energy in the
UV wavelength range. In Fig. 2, we compare also the NLTE derivatives with the
LTE ones. The largest difference can be seen again in the $U$ passband.

Effect of NLTE on the flux derivatives as a function of effective temperature
for the Strömgren $uy$ passbands is presented in Fig. 3. In general,
differences between the LTE and NLTE derivatives are relatively small and they
increase with effective temperature, especially in the case of $\alpha_{g}$ in
the $u$ passband. In the case of both model atmospheres, the values of
$\alpha_{T}$ in the $u$ passband change only slightly with the temperature,
whereas in the $y$ passband they are larger for higher $T_{\rm eff}$. The
gravity derivatives get more steep with the higher temperature, except
$\alpha_{g}^{u}(T_{\rm eff})$ for NLTE models, which is a decreasing function
up to $T_{\rm eff}\approx 24000$ K and then it increases.

Figure 1: The NLTE flux derivatives over $\log T_{\rm eff}$ (the left panel)
and $\log g$ (the right panel) as a function of temperature for the Strömgren
$uy$ passbands. Three values of $\log g=3.5,4.0,4.5$ were considered. The
values of the atmospheric metallicity and microturbulent velocity were assumed
as $(Z/Z_{\odot})_{\rm atm}$=1 and $\xi_{t}$=2 km/s, respectively. Figure 2: A
comparison of NLTE and LTE flux derivatives as a function of the central
wavelengths of the $UBVRIJHK$ passbands, for a models with $T_{\rm eff}$=20000
K, $\log g=4.0$, $(Z/Z_{\odot})_{\rm atm}$=1, $\xi_{t}$=2 km/s. The left and
right panel correspond to the temperature ($\alpha_{T}^{x}$) and gravity
($\alpha_{g}^{x}$) derivatives, respectively. Figure 3: Values of
$\alpha_{T}^{x}$ (the left panel) and $\alpha_{g}^{x}$ (the right panel) as a
function of temperature for LTE and NLTE models and the Strömgren $uy$
passbands. The remaining parameters of the model are: $\log g=4.0$,
$(Z/Z_{\odot})_{\rm atm}$=1, $\xi_{t}$=2 km/s.

###  3.2 Limb darkening

Knowledge of a distribution of the specific intensity over the stellar disc is
crucial in many fields of astrophysics. In the case of pulsating stars, limb-
darkening function and its derivatives over effective temperature and gravity,
are the second atmospheric data needed to compute theoretical values of the
photometric amplitudes and phases.

A usage of analytical formulae for the specific intensity distribution saves
enormously the computation time. The first limb-darkening law was proposed by
Milne (milne (1921)) in a linear form. With the development of models of
stellar atmospheres, it turned out that this approximation has a poor
accuracy. Thereafter, Klinglesmith & Sobieski (kinglesmith (1970)) published a
logarithmic law for early type stars. A quadratic formula for the limb-
darkening law was proposed by Manduca at al. (manduca (1977)) and Wade &
Ruciński (1985). For stars hotter than 8500 K, Diáz-Cordovés & Giménez (diaz
(1992)) suggested a square root law.

However, all these above formulas were not adequate for all types of stars. A
more general law was proposed by Claret (claret2000 (2000)) in the following
non-linear form

$\frac{I(\lambda,\mu)}{I(\lambda,1)}=1-\sum_{k=1}^{4}a^{\lambda}_{k}(1-\mu^{\frac{k}{2}})$
(8)

where $I(\lambda,\mu)$ is a specific intensity at the wavelength $\lambda$,
$\mu$ is a cosine of the angle between a line of sight and the emergent
intensity, and $I(\lambda,1)$ is the value at the center of the stellar disc.
Limb darkening coefficients, $a^{\lambda}_{k}$, are determined to reproduce
the model intensity distribution and to conserve the flux with a high accuracy
in the whole range of effective temperatures and gravities.

These coefficients were calculated for LTE model atmospheres in many
photometric systems and a wide range of effective temperatures, gravities,
metallicities and for several values of the microturbulent velocities (Claret
2000, 2003, Claret & Hauschildt 2003, Claret 2008). A comparison of different
numerical methods of computations of limb darkening coefficients (LDC) were
discussed in detail by Díaz-Cordovés et al. (diaz1995 (1995)), Claret
(claret2000 (2000)), Heyrovský (heyrovsky (2007)) and Claret (claret2008
(2008)).

In this paper, we determine non-linear LDC for metal line-blanketed, NLTE,
plane-parallel, hydrostatic model atmospheres of Lanz & Hubeny (lanzhubeny
(2007)). In the first step, we computed monochromatic specific intensities for
20 equally separated points of $\mu$ in the range of $\left<0.001,1\right>$
and in the wavelength range of (2950,26500) $\AA$, for all parameters of the
BSTAR2006 atmospheres. All calculations were performed using the SYNSPEC
program (Hubeny at al. synspec2 (1985), Hubeny & Lanz synspec (2000)). This
code is intended to compute specific intensities and fluxes from the model
atmosphere input, which we took from NLTE models described above.
Subsequently, for each disc position, $\mu$, we integrated intensities over
the $uvby$ and $UBVRIJHK$ bands

$I_{x}(\mu)=\frac{\int\limits^{\lambda_{1}}_{\lambda_{2}}I(\lambda,\mu)S(\lambda)\,d\lambda}{\int\limits^{\lambda_{1}}_{\lambda_{2}}S(\lambda)\,d\lambda},$
(9)

where $I_{x}(\mu)$ is the specific intensity in the passband $x$ and
$I(\lambda,\mu)$ is the monochromatic specific intensity obtained by SYNSPEC.
The response function, $S(\lambda)$, was interpolated to the wavelengths of
computed monochromatic specific intensity by the cubic spline method.

The non-linear limb darkening coefficients, $a_{k}^{x}$, were determined using
the least squares method by minimizing

$\chi^{2}(a_{1},a_{2},a_{3},a_{4})=\sum_{i=1}^{20}\left(I_{LDC_{i}}-I_{i}\right)^{2},$
(10)

where $I_{LDC_{i}}$ is the specific intensity computed with the Claret
nonlinear law and $I_{i}$ is the corresponding model intensity, at the point
$\mu_{i}$.

Similarly to Heyrovský (heyrovsky (2007)), we evaluated the quality of our
fits by the relative residuals

$\sigma^{2}=\frac{\sum\limits_{i=1}^{20}\left(I_{i}-I_{LDC_{i}}\right)^{2}}{\sum\limits_{i=1}^{N}I_{i}^{2}},$
(11)

whereas, the conservation of the flux was controlled by computing the relative
flux excess

$\frac{\Delta{\cal F}_{x}}{{\cal
F}_{x}}=\frac{\int\limits_{0}^{1}I_{x}\mu\,d\mu-\int\limits_{0}^{1}I_{x,LDC}\mu\,d\mu}{\int\limits_{0}^{1}I_{x}\mu\,d\mu}.$
(12)

The limb darkening law, $h(\mu)$, given in Eq. 6, is defined as (Daszyńska-
Daszkiewicz et al. 2002)

$h_{x}(\mu)=2\pi\leavevmode\nobreak\ \frac{I_{x}(\mu)}{{\cal
F}_{x}}=2\pi\leavevmode\nobreak\
\frac{1-\sum\limits_{k=1}^{4}a_{k}^{x}(1-\mu^{k/2})}{1-\sum\limits_{k=1}^{4}\frac{k}{k+4}a_{k}^{x}}.$
(13)

Names of the LDC files are given Table 1 and they are coded in the same way as
the flux and flux derivatives tables (cf. Sect. 3.1). Each file contains the
following columns: names of LDC, $a_{k}$, effective temperature, $T_{\rm
eff}$, gravity, $\log g$, metallicity, $(Z/Z_{\odot})_{\rm atm}$,
mictroturbulent velocity, $\xi_{t}$, and values of LDC in the Strömgren and
Johnson-Cousins-Glass photometric passbands, in the order: $uvbyUBVRIJHK$.

In Fig. 4 we show distributions of the normalized specific intensity in the
$UBVRIJHK$ passbands for a model with the following parameters: $T_{\rm
eff}$=20000 K, log g=4.0, $\xi_{t}$=2 km/s, $(Z/Z_{\odot})_{\rm atm}$=1. In
this figure, we compare the specific intensity computed by means of SYNSPEC
for 20 equally separated points of $\mu$ (squares) and the fitted limb
darkening function defined by Eq. 8 (solid lines).

Figure 4: The angular distribution of the normalized specific intensity for
the $UBVRIJHK$ passbands and the NLTE model with $T_{\rm eff}$=20000 K, $\log
g$=4.0, $\xi_{t}$=2 km/s and $(Z/Z_{\odot})_{\rm atm}$=1. The actual model
intensities are marked as squares and intensities computed from LDC as solid
lines. For clarity, passbands from $B$ are shifted downward by 0.1.

Intensities for passbands $BVRIJHK$ were shifted downward by $n\cdot 0.1$,
where n=1,2,…,7. As we can see a quality of the limb darkening fit is very
accurate.

To compare values of the NLTE and LTE intensities, in Fig. 5 we plot $I_{x}$
[erg$\cdot$cm${}^{-2}\cdot$s${}^{-1}\cdot$ster-1] as a function of the angle
$\mu$ for the Strömgren $uvby$ passbands. The same model as in Fig. 4 was
considered. The NLTE intensities get lower values, in particular, for the $u$
passband. This is caused by the location of the $u$ filter on the Balmer
continuum where the difference in the amount of energy radiated in LTE and
NLTE models is more pronounced than in a region to the right from the Balmer
jump, where other filters ($vby$) are defined.

Figure 5: The intensity distribution for the $uvby$ passbands for the NLTE
(solid lines) and LTE (dashed lines) models with $T_{\rm eff}$=20000 K, $\log
g$=4.0, $\xi_{t}$=2 km/s and $(Z/Z_{\odot})_{\rm atm}$=1.

In Fig. 6 and 7, we plot the NLTE limb darkening coefficients,
$a_{1},a_{2},a_{3},a_{4}$, as a function of $T_{\rm eff}$ in the Strömgren $u$
passband, for different values of gravity, $\log g$, and metallicity,
$(Z/Z_{\odot})_{\rm atm}$, respectively. All panels have the same scale. As we
can see, there is a strong dependence of LDC on effective temperature, gravity
and metallicity. The sensitivity on $T_{\rm eff}$ is stronger for lower values
of $\log g=3.5$. The sensitivity to metallicity is similar for all LDC. The
effect of the microturbulent velocity, $\xi_{t}$, not shown here, is
comparable to the effect of metallicity.

Figure 6: The NLTE limb darkening coefficients in the $u$ passband as a
function of $T_{\rm eff}$ for three values of gravity: $\log=3.5$ (solid
line), $\log=4.0$ (dashed line) and $\log=4.5$ (dotted line). The solar
metallicity, $(Z/Z_{\odot})_{\rm atm}$=1, and microturbulent velocity of
$\xi_{t}$=2 km/s were assumed. Figure 7: The NLTE limb darkening coefficients
in the $u$ passband as a function of $T_{\rm eff}$ for six values of the
atmospheric metallicity: $(Z/Z_{\odot})_{\rm atm}$=2.0 (solid black line), 1.0
(dashed black line), 0.5 (dotted black line), 0.2 (solid grey line), 0.1
(dashed grey line) and 0.0 (dotted grey line). The gravity of $\log g$=4.0 and
microturbulent velocity of $\xi_{t}$=2 km/s were assumed.

Let us now discuss the accuracy of our determinations of the Claret limb-
darkening coefficients for the BSTAR2006 NLTE models. In Table. 2, we give the
average and maximum values of the merit functions, i.e. the relative
residuals, $\sigma$, and relative flux excess, $\left|\Delta F\over F\right|$.
These values were calculated for the whole range of $T_{\rm eff}$, i.e.
(15000,30000) K, and $\log g$, i.e. (1.75, 4.75) dex, assuming the solar
metallicity, $(Z/Z_{\odot})_{\rm atm}$=1, and the microturbulent velocity of
$\xi_{t}$=2 km/s. We used 20 equally spaced points of $\mu$ instead of 17
unequally spaced points as used in LTE Kurucz models. Moreover, our NLTE
values of $I(\mu)$ were computed closer to the stellar limb, i.e. up to
$\mu=0.001$, whereas the lowest angle of $\mu$ in LTE models is 0.01. We
considered two sets of passbands: $uvbyUBVRIJHK$ and $BVRI$. The second set
was used to compare our results with those of Heyrovský (heyrovsky (2007)) who
made computations for LTE Kurucz atmosphere models using 17 and 11 points of
$\mu$. As we can see the quality of our fit is very good. The average value of
$\sigma$ amounts to 0.119% and for the worst fitted profile we got
$\sigma$=0.432%. Our average value of $\sigma$ is smaller than the Heyrovský’s
one by a factor of two. When we limited our analysis to the $BVRI$ passbans,
as in Heyrovský (heyrovsky (2007)), the result is even better. In this case,
the average and maximum values of $\sigma$ amount to 0.089% and 0.251%,
respectively. This indicates that with the larger number of equally separated
points of $\mu$ in the fitting procedure one reproduces more accurately the
model intensities.

Table 2: Comparison of our limb-darkening fit quality with results of
Heyrovsky (2007) for the Claret non-linear law. In columns 2 and 3, we give
the average and maximum values of the relative rms residual, $\sigma$,
respectively. Columns 4 and 5 contain the average and maximum values of the
relative flux excess, $\left|\Delta F\over F\right|$. Passbands for which
these goodness-of-fit quantities were evaluated are given in notes.

Method | Average $\sigma$ | Max. $\sigma$ | Average $\left|\Delta F\over F\right|$ | Max. $\left|\Delta F\over F\right|$
---|---|---|---|---
| % | % |  |
20-pointa | 0.119 | 0.432 | $14.5\times 10^{-5}$ | $9.10\times 10^{-4}$
20-pointb | 0.089 | 0.251 | $9.50\times 10^{-5}$ | $5.28\times 10^{-4}$
17-pointc | 0.190 | 0.625 | $8.27\times 10^{-5}$ | $5.05\times 10^{-4}$
11-pointc | 0.175 | 0.558 | $4.57\times 10^{-5}$ | $2.17\times 10^{-4}$
222 a values evaluated for all 12 passbands: $uvbyUBVRIJHK$. b values
evaluated for the $BVRI$ passbands. c values from Heyrovský (heyrovsky (2007))
for the $BVRI$ passbands.

Finally, let us check the conservation of the flux. The average value of
$\left|\Delta F\over F\right|$ from our fitting for NLTE models is slightly
worse than in Heyrovský (heyrovsky (2007)) for LTE models atmospheres, but of
the same order of magnitude. The maximum values of $\left|\Delta F\over
F\right|$ are $9.10\times 10^{-4}$ for all 12 passbnds and $5.28\times
10^{-4}$ for $BVRI$. The corresponding average values are $14.5\times 10^{-5}$
and $9.50\times 10^{-5}$, respectively. For larger values of gravities, the
fluxes computed with the Claret limb darkening law are overestimated whereas
for smaller values of $\log g$, fluxes are underestimated. This is caused by a
steeper slope of the intensity near the limb for lower values of $\log g$.
Consequently, the fitted limb darkening for lower gravities is slightly below
model intensities.

##  4 Uncertainties in photometric pulsational observables from model
atmospheres

Figure 8: A comparison of the NLTE and LTE values of the amplitude ratios,
$A_{u}/A_{y}$ (the left panel) and phase differences,
$\varphi_{u}-\varphi_{y}$ (the right panel) for 10 $M_{\odot}$ main sequence
models as a function of $T_{\rm eff}$ for the first three radial modes: p1,
p2, p3. The atmospheric metallicity of $(Z/Z_{\odot})_{\rm atm}=1$ and
microturbulent velocity of $\xi_{t}=2$ km/s were assumed. Linear nonadiabatic
pulsations were computed with hydrogen and metal abundance of $X=0.7$ and
$Z=0.02$, respectively, the OPAL opacities and the AGSS09 chemical mixture.
Figure 9: The same as in Fig. 8 but for two $\ell=1$ modes: p1 and g14. Figure
10: The same as in Fig. 8 for for two $\ell=2$ modes: p1 and g22. Figure 11:
Effect of the atmospheric metallicity on the amplitude ratios, $A_{u}/A_{y}$
(the left panel) and phase differences, $\varphi_{u}-\varphi_{y}$ (the right
panel) for 10 $M_{\odot}$ main sequence models as a function of $T_{\rm eff}$
for three radial modes: p1, p2, p3. NLTE model atmospheres with the
microturbulent velocity of $\xi_{t}=2$ km/s and two values of
$(Z/Z_{\odot})_{\rm atm}$ were used.

The atmospheric metallicity, $(Z/Z_{\odot})_{\rm atm}$, microturbulent
velocity, $\xi_{t}$, and effects of NLTE are the most important factors which
can affect the photometric amplitudes and phases of a pulsating star.

In all comparisons, we used a reference model computed with the mass of
$M=10M_{\odot}$, the OPAL opacity tables, the AGSS09 mixture, hydrogen
abundance of $X$=0.7, metal abundance of $Z$=0.02, without overshooting from a
convective core, $\alpha_{\rm ov}$=0.0, and NLTE-TLUSTY models of stellar
atmospheres with the metallicity of $(Z/Z_{\odot})_{\rm atm}$=1.0 and the
microturbulent velocity of $\xi_{t}$=2 km/s.

In Fig. 8, we compare the photometric observables computed with LTE and NLTE
model atmospheres at the same values of metallicity, $(Z/Z_{\odot})_{\rm
atm}=1$, and microturbulent velocity, $\xi_{t}=2$ km/s, for the $10M_{\odot}$
model in the course of its main sequence evolution. We considered the
Strömgren $uy$ passbands and the first three radial modes: p1, p2, p3. In the
left panel, we show the amplitude ratios, $A_{u}/A_{y}$, and in the right one
the corresponding phase differences, $\varphi_{u}-\varphi_{y}$, as a function
of $T_{\rm eff}$. As we can see, the fundamental mode is most sensitive to the
NLTE effects. In general, the NLTE values of the amplitude ratio,
$A_{u}/A_{y}$, are smaller than the LTE ones for the hotter models and larger
for more evolved models.

Subsequently we checked the NLTE effect for nonradial modes. In Fig. 9, we
plot the same photometric observables as in Fig. 8, but for the $\ell=1$ mode
considering one pressure mode, p1, and one high-order gravity mode, g14. Fig.
10 illustrates the same but for two $\ell=2$ modes: p1 and g22. In the case of
nonradial modes, the amplitude ratios, $A_{u}/A_{y}$, have larger values for
NLTE computations. We can see also a different behavior and values of
$A_{u}/A_{y}$ and $\varphi_{u}-\varphi_{y}$ for pressure and high-order
gravity modes, especially for the $\ell=2$ modes.

Effect of the atmospheric metallicity on the photometric observables for the
radial modes are shown in Fig. 11. Here, we considered $(Z/Z_{\odot})_{\rm
atm}=0.1$ in addition to our standard value of $(Z/Z_{\odot})_{\rm atm}=1.0$.
As we can see, the effect of the atmospheric metallicity is comparable to the
effect of the departure from LTE in stellar atmospheres. In the case of
nonradial modes the effect of $(Z/Z_{\odot})_{\rm atm}$ is negligible.

Influence of the microturbulent velocity, $\xi_{t}$, is of the same order as
$(Z/Z_{\odot})_{\rm atm}$, as has been shown by Szewczuk & Daszyńska-
Daszkiewicz (2010).

##  5 Uncertainties in pulsational photometric observables from pulsation
theory

Linear nonadiabatic theory of stellar pulsation provides eigenfrequencies,
corresponding eigenfunctions and information on mode excitation. The value of
the flux eigenfunction at the level of the photosphere enters the expression
for the light variation and it is called the $f$-parameter, as it has been
introduced in Sect. 2. For a given mode frequency, this parameter depends, in
general, on mean stellar parameters, chemical composition, opacity data,
subphotospheric convection etc. Moreover, in the case of high-order gravity
modes, the $f$-parameter depends also on the mode degree, $\ell$ (e.g.
Daszyńska-Daszkiewicz & Walczak 2010). There are also large differences in
values of $f$ between high and low frequency modes of a given degree, $\ell$.
This is why the amplitude ratios and phase differences for p- and high-order
g-modes behave so different, which is exemplified most pronounced in Fig. 10.

Figure 12: Effects of the internal abundances of hydrogen, $X$, and metals,
$Z$, on photometric observables for the first three radial modes as a function
of $T_{\rm eff}$ for the 10 $M_{\odot}$ main sequence models. The NLTE models
with the atmospheric metallicity of $(Z/Z_{\odot})_{\rm atm}=1$ and
microturbulent velocity of $\xi_{t}=2$ km/s were used. Linear nonadiabatic
pulsations were computed with OPAL opacity tables and AGSS09 chemical mixture.
Figure 13: The same as in Fig. 12 but for two $\ell=1$ modes: p1 and g14.
Figure 14: The same as in Fig. 12 but for two $\ell=2$ modes: p1 and g22.
Figure 15: Effects of the opacity tables and chemical mixture on photometric
observables for the first three radial modes as a function of $T_{\rm eff}$
for the 10 $M_{\odot}$ main sequence models. The NLTE models with the
metallicity of $(Z/Z_{\odot})_{\rm atm}=1$ and microturbulent velocity of
$\xi_{t}=2$ km/s were used. Linear nonadiabatic pulsations were computed with
$X=0.7$ and $Z=0.02$. Figure 16: The same as in Fig. 15 but for two $\ell=1$
modes: p1 and g14. Figure 17: The same as in Fig. 15 but for two $\ell=2$
modes: p1 and g22.

Pulsations of the B-type stars are driven by the $\kappa$ mechanism operating
in the ,,Z-bump” layer. The nonadiabatic complex parameter, $f$, is also
defined in this subphotospheric layer, therefore its value strongly depends on
metal abundances, chemical mixture and hence opacities (Daszyńska-Daszkiewicz,
Dziembowski, Pamyatnykh 2005).

In Fig. 12, we show effects of the internal abundances of hydrogen, $X$, and
metals, $Z$, on the photometric amplitude ratios, $A_{u}/A_{y}$, and phase
differences, $\varphi_{u}-\varphi_{y}$, for the first three radial modes: p1,
p2, p3. The same model as in Sect. 4 was selected. Here, we compare
computations obtained with $Z$=0.02 $vs.$ $Z$=0.03 and $X$=0.7 $vs.$ $X$=0.75.
As we can see, the amplitude ratios computed with $Z$=0.03 are larger than
those obtained with $Z$=0.02, whereas increasing hydrogen abundance, $X$,
decreases the amplitude ratios. The value of $\varphi_{u}-\varphi_{y}$ changes
by about 0.1 rad at most. These effects are particulary distinct for the
radial fundamental mode, p1. Comparing Fig. 12 with Fig. 8 and 11, we can see
also that the effects of $X$ and $Z$ on the photometric observables are much
more pronounced than the effects of NLTE and $(Z/Z_{\odot})_{\rm atm}$. Please
note that for a given mode degree, $\ell$, scales of the amplitude ratios and
phases differences are the same in all figures.

Effects of the hydrogen and metallicity abundance for the $\ell=1$, p1 and
$g_{14}$ modes is shown in Fig. 13, and for the $\ell=2$, p1 and g22 modes in
Fig. 14. As one can see, in the case of the nonradial modes these parameters
have smaller influence on the photometric observables than in the case of
radial modes. The amplitude ratio, $A_{u}/A_{y}$, for the $\ell=1$ mode
increases with increasing abundances of both $Z$ and $X$. This is true also
for the $\ell=2$, p1 mode, whereas the amplitude ratio of the $\ell=2$, g22
mode does not behave monotonically. The important result is that in the case
of nonradial modes the effects of $X$ and $Z$ are considerable smaller than
the NLTE effects (cf. Fig. 9 and 10), which is opposite to the radial mode
case.

Finally, we evaluated effects of the opacity data and chemical mixture. In
Fig. 15 we show the photometric observables for the radial modes computed with
two sources of the opacity tables and two chemical mixtures. Here, we compared
computations obtained with the OPAL $vs.$ OP tables and the GN93 $vs.$ AGSS09
mixture. Again, the largest effects are for the p1 mode in both the amplitude
ratios and phase differences. The amplitude ratios computed with GN93 are
smaller than those obtained with AGSS09, because of relatively higher
abundance of iron in AGSS09. Also computations with the OP data give smaller
values of $A_{u}/A_{y}$ comparing to those obtained with the OPAL data.
Moreover, the maximum of the OP amplitude ratio is shifted towards higher
effective temperature. This is connected with the location of the $Z-$bump
layer, which occurs at slightly higher temperature in the OP data. A usage of
different opacity data changes the value of $\varphi_{u}-\varphi_{y}$ by about
0.2 rad at most. The chemical mixture has rather a minor effect on phase
differences. We can see that also these input data affect the photometric
observables of radial modes far more than the atmospheric parameters (cf. Fig.
8 and 11).

Effects of the opacity tables and chemical mixture for nonradial modes with
$\ell=1$ and 2 are shown in Fig. 16 and Fig. 17, respectively. These data
affect the $\ell=1$ modes very subtle and are far less important than the NLTE
effects (cf. Fig. 9). In the case of the $\ell=2$ modes effects of opacities
and mixture are comparable to the NLTE effects (cf. Fig. 10).

As we could see, generally, in the case of radial modes effects of pulsational
parameters on photometric observables are much larger than effects of
atmospheric parameters discussed in previous section. The opposite is true for
nonradial modes which are more sensitive to the NLTE effects, in particular in
the case of dipole modes.

##  6 Conclusions

We have presented a comprehensive overview of possible sources of
uncertainties in theoretical values of the photometric amplitudes and phases
of early B-type pulsators. These data are of particular importance because
they serve as tools for mode identification in the case of main sequence
pulsators, as well as contain information on stellar physics.

The uncertainties are embedded in stellar model atmospheres and nonadiabatic
theory of stellar pulsation. The atmospheric input consists of the flux
derivatives over $T_{\rm eff}$ and $\log g$, and limb darkening and its
derivatives, in photometric passbands. These quantities are sensitive to the
atmospheric metallicity and microturbulent velocity, as well as to the
departure from the LTE approximation. From pulsation computations, we get the
nonadiabatic complex parameter, $f$, whose value depends on chemical
composition, opacities and subphotospheric convection if present. If the
$f$-parameter can be derived from observations then we get an extra seismic
probe by means of which we can test this input physics.

We begun with computations of tables with various quantities, needed to
evaluate the light variation in a given photometric band, for NLTE-TLUSTY
model atmospheres. These tables include data on the passband flux, flux
derivatives over effective temperature and gravity, and the non-linear limb
darkening coefficients in 12 most popular passbands: $uvbyUBVRIJHK$. All these
data are public available and can be retrieved from the Wrocław HELAS web
page.

Then, we studied effects of these input parameters on the photometric
observables, i.e. amplitude ratios and phase differences. We considered the 10
$M_{\odot}$ main sequence models and the low degree modes with $\ell=0,1,2$.
In particular, the effect of NLTE model atmospheres was studied for the first
time. In the case of radial modes, this effect is comparable to an influence
of the atmospheric metalicity and microturbulent velocity. The photometric
observables of nonradial modes appeared more sensitive to the departure from
LTE. Subsequently, we drew a parallel to effects of parameters related to
pulsation computations. These comparisons showed that the photometric
observables of the radial modes are by far more sensitive to the pulsation
input. In turn, in the case of nonradial modes the NLTE effect became more
important, especially for the $\ell=1$ modes.

A complete seismic model should reproduce not only pulsational frequencies but
also the observed values of the photometric amplitude and phases, which can be
translated into the empirical values of the $f$-parameter. Our studies showed
that such comprehensive approach should take into account all inaccuracies in
the photometric amplitudes and phases.

###### Acknowledgements.

This work was supported by the HELAS EU Network, FP6, No. 026138 and for JDD
by the Polish MNiSW grant N N203 379636.

## References

  * (1) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481 (AGSS09)
  * (2) Balona L. A., Stobie R. S. 1979, MNRAS, 189, 649
  * (3) Claret A. 2000, A&A, 363, 1081
  * (4) Claret A. 2003, A&A, 401, 657
  * (5) Claret A. 2008, A&A, 482, 259
  * (6) Claret A., Hauschildt P. H. 2003, A&A, 412, 241
  * (7) Cugier H., Dziembowski W. A, Pamyatnykh A. A. 1994, A&A, 291, 143
  * (8) Daszyńska-Daszkiewicz J., Dziembowski W. A., Pamyatnykh A. A., Goupil M. J. 2002, A&A, 392, 151
  * (9) Daszyńska-Daszkiewicz J., Dziembowski W. A., Pamyatnykh A. A. 2003, A&A, 407, 999
  * (10) Daszyńska-Daszkiewicz J., Dziembowski W. A., Pamyatnykh A. A. 2005, A&A, 441, 641
  * (11) Daszyńska-Daszkiewicz J., Dziembowski W. A., Pamyatnykh A. A. 2007, AcA, 57, 11
  * (12) Daszyńska-Daszkiewicz J., Walczak P., 2010, MNRAS, 403, 496
  * (13) Diáz-Cordovés J., Giménez A. 1992, A&A, 259, 227
  * (14) Diáz-Cordovés J., Claret A., Giménez A. 1995, A&AS, 110, 329
  * (15) Dziembowski W. A. 1977, AcA, 27, 203
  * (16) Grevesse N., Sauval A. 1998, Space Sci. Rev., 85, 161 (GS98)
  * (17) Grevesse N., Noels A. 1993, in Origin ans Evolution of the Elements, ed. N. Pratzo, E. Vangioni-Flam, & M. Casse (Cambridge Univ. Press), 15 (GN93)
  * (18) Heyrovský D. 2007, AJ, 656, 492
  * (19) Hubeny I., Lanz T., 1995, ApJ, 439, 875
  * (20) Hubeny I., Lanz T. 2000, SYNSPEC - A User’s Guide
  * (21) Hubeny I., Stefl S., Harmanec P. 1985, Bull. Astron. Inst. Czechosl., 36, 214
  * (22) Iglesias C. A., Rogers F. J. 1996, ApJ, 464, 943
  * (23) Kinglesmith D. A., Sobieski S. 1970, AJ, 75, 175
  * (24) Kurucz, R. L. 2004, http:// kurucz.harvard.edu
  * (25) Lanz T., Hubeny I. 2003, ApJS, 146, 417
  * (26) Lanz T., Hubeny I. 2007, ApJS, 169, 83
  * (27) Manduca A., Bell R. A., Gustafsson B. 1977, A&A, 61, 809
  * (28) Milne E. A. 1921, MNRAS, 81, 361
  * (29) Moro D., Munari U. 2000, A&AS, 147, 361
  * (30) Seaton M. J. 2005, MNRAS, 362, L1
  * (31) Stamford P. A., Watson R. D. 1981, Ap&SS, 77, 131
  * (32) Szewczuk W., Daszyńska-Daszkiewicz J. 2010, AN, in press
  * (33) Townsend R. 2003, MNRAS, 343, 125
  * (34) Wade R. A., Ruciński S. M. 1985, A&AS, 60, 471
  * (35) Watson R. D. 1988, Ap&SS, 140, 255