# Spin-$1/2$ invisible particles in heavy meson decays

Geng<EMAIL_ADDRESS>Tianhong<EMAIL_ADDRESS>(Corresponding author), Yue<EMAIL_ADDRESS>Jing-Bo
<EMAIL_ADDRESS>and Guo-Li<EMAIL_ADDRESS>
1School of Physics, Harbin Institute of Technology, Harbin, 150001, China
2Department of Physics, Hebei University, Baoding 071002, China
3Hebei Key Laboratory of High-Precision Computation and Application of Quantum
Field Theory, Baoding 071002, China

###### Abstract

The flavor-changing neutral current decay processes of the $B$ and $B_{c}$
mesons with the final states involving spin-$1/2$ particles are investigated.
By considering the background of the Standard Model where $\nu\bar{\nu}$
contributes the missing energy and the experimental upper bounds for the
branching fractions, we get the constraints of the coupling constants of the
quark-antiquark and the assumed invisible particles $\chi\bar{\chi}$. The
constraints of the coupling constants are then used to study the similar
processes of the $B_{c}$ meson. At some specific region of $m_{\chi}$, the
upper limit of BR($B_{c}\to D_{(s)}\chi\bar{\chi}$) is of the order of
$10^{-6}$, while for BR($B_{c}\to D^{\ast}_{(s)}\chi\bar{\chi}$), it is
$10^{-5}$. The possibility of distinguishing $\chi$ to be a Majorana or Dirac
fermion by the differential branching fractions is also discussed.

##  I Introduction

As the freeze-out mechanism Bernstein  _et al._ (1985); Srednicki  _et al._
(1988) can naturally interpret the observed dark matter abundance in our
Universe, the weakly interacting massive particle (WIMP) is considered to be
one of the most promising dark matter candidates. It is considered as a
thermal relic from the local thermodynamic equilibrium early Universe
Izaguirre  _et al._ (2015). The observed dark matter relic abundance
$\Omega_{c}h^{2}=0.1131\pm 0.0034$ Bertone  _et al._ (2005); Komatsu  _et al._
(2009) sets a lower bound for WIMP’s annihilation cross section. In specific
models, the cross section can be connected to the mass of WIMP and coupling
coefficients between WIMP and the Standard Model (SM) fermions. For example,
the Lee-Winberg limit Lee and Weinberg (1977) demands its mass larger than a
few GeV. However, this result is model dependent. With different models or
proper selection of parameters, this constraint can be relaxed, which makes
lower mass WIMP to be possible. For example, the MeV-scale light dark matter
(LDM) was proposed Pospelov  _et al._ (2008); Hooper and Zurek (2008) to
explain the unexpected emission of 511 keV photons from the galaxy center.

Previous experiments mostly focused on the dark matter particle with large
mass, namely hundreds of GeV to several TeV. But recent experiment Aprile  _et
al._ (2017) has set much stricter constraints on the parameter space for the
WIMP with mass larger than several GeV. It provides a motivation to study the
sub-GeV LDM through high-energy colliders. For example, CODEX-b at the LHCb
experiment aimed to probe for GeV-scale long-lived particles Gligorov  _et
al._ (2018). Missing energy signals Abdallah  _et al._ (2014) in flavor-
changing neutral current (FCNC) processes of heavy mesons provide a possible
way to probe light WIMP. Within the SM, neutrinos $\nu\bar{\nu}$ in the final
state make contribution to the missing energy. However, theoretical
calculations of the branching fractions of $B\to h_{f}\nu\bar{\nu}$ are less
than the experimental bounds of $B\to h_{f}\not{E}$, where $h_{f}$ is the
final meson and $\not{E}$ is the missing energy. So there is still some
allowed parameter space for the decays involving other light invisible
particles.

Theoretically, spin of the invisible particle has several possibilities
Belyaev  _et al._ (2019). It can be a (pseudo)scalar Boehm and Fayet (2004), a
fermion Kusenko (2009), or a hidden vector Hambye (2009). In the previous
paper Li  _et al._ (2019), we have considered the scalar and pseudoscalar
cases. In this paper, we focus on the spin-$1/2$ light dark matter particles.
There are many models involving the fermionic dark matter particles, such as
sterile neutrino Kusenko (2009), neutralino Drees and Hajkarim (2018); Yue and
Chu (2018), Higgs-portal Djouadi  _et al._ (2012, 2013), Z-portal Arcadi  _et
al._ (2015) and singlet-doublet Hisano  _et al._ (2018); Barman  _et al._
(2019); Restrepo  _et al._ (2019); Abe and Sato (2019); Fiaschi  _et al._
(2019). Specifically, it can be either a Majorana or a Dirac fermion, as it is
electrically neutral. Phenomenologically, the new invisible fermion can weakly
interact with the SM fermions via a mediator, which can be a scalar Matsumoto
_et al._ (2019), pseudoscalar Yang (2016), vector or axial-vector Chala  _et
al._ (2015) particle. The mass of the mediator is usually considered to be
hundreds of GeV. In the energy level of heavy meson decays, namely several
GeV, the branching ratios are greatly suppressed. However, as the FCNC and
annihilation processes in the SM are also highly suppressed, the contribution
of the new physics maybe important, which has been extensively studied in the
decays of mesons Bird  _et al._ (2004, 2006); Badin and Petrov (2010);
Gninenko and Krasnikov (2015); Barducci  _et al._ (2018); Kamenik and Smith
(2012); Bertuzzo  _et al._ (2017). For example, Ref. Badin and Petrov (2010)
most focused on the $B$ meson annihilation processes, and Ref. Kamenik and
Smith (2012) studied various dark sectors in $B$ meson FCNC processes.

In this work, we will further study the spin-1/2 invisible particles in $B$
and $B_{c}$ meson FCNC decays. Such studies for the $B_{c}$ meson are still
missing. The $B_{c}$ meson has been massively produced and measured by the CDF
Aaltonen  _et al._ (2016), ATLAS Burdin (2016), CMS Berezhnoy  _et al._
(2019), and LHCb Aaij  _et al._ (2019) experiments. The production rate of
$B_{c}$ meson on the LHCb Collaboration is close to 3.7 per mille of that of
the $B$ mesons Aaij  _et al._ (2019). The $B_{c}$ events are in the order of
$10^{10}$ per year. As the luminosity of the LHC increases significantly, much
more $B_{c}$ events will be generated in the near future. We will first
introduce the effective operators to describe the coupling between quarks and
the invisible fermions. The experimental upper bounds for the FCNC decay
channels of the $B$ meson will then provide constraints of the coupling
constants, which will be applied to calculate the upper bounds of the similar
decay processes of the $B_{c}$ meson. To calculate the hadronic transition
matrix elements, two methods are used: for the $B\rightarrow h_{f}$, the QCD
light-cone sum rules (LCSR) is used, while for the $B_{c}\rightarrow h_{f}$,
we apply the instantaneous Bethe-Salpeter (BS) method which is more suitable
for such cases. For the light invisible fermions, both the Majorana and Dirac
cases are considered. As they interact differently with quarks, for example,
the Majorana fermion has neither vector nor tensor interactions, while Dirac
fermion has both of them; the differential distribution will show slight
difference.

The paper is organized as follows: in Sec. II, we present the model-
independent effective Lagrangian to describe the coupling between the light
invisible fermions and quarks, and extract the constraints of the coupling
coefficients. In Sec. III, we calculate the upper limits of the branching
fractions of $B_{c}$ decays, and give the differential decay rate as a
function of the missing energy. Finally, we draw the conclusion in Sec. IV.

##  II Effective operators

The FCNC decay processes of heavy meson to spin-$1/2$ invisible particles
$\bar{\chi}\chi$ are described in Fig. 1, where $q$, $q_{{}_{f}}$, and
$\bar{q}^{\prime}$ represent the quark and antiquark, respectively.

Figure 1: Feynman diagrams of decay channels involving invisible particles.

The four-fermion vertex may be generated at the tree or loop level by
introducing new physical mediators in specific models Matsumoto  _et al._
(2019); Yang (2016); Chala  _et al._ (2015). In this work, we follow Ref.
Badin and Petrov (2010) to introduce a model-independent effective Lagrangian,

$\displaystyle\mathcal{L}_{eff}=\sum_{i=1}^{9}{g_{{}_{fi}}Q_{i}},$ (1)

where the subscript $f$ represents fermion which can be Majorana or Dirac
type, and $g_{fi}$ are the phenomenological coupling constants which are
suppressed by the new physical energy scale $\Lambda$. There are nine
independent effective operators $Q_{i}$s, which have the forms

$\displaystyle
Q_{1}=(\bar{q}_{{}_{f}}q)(\bar{\chi}\chi),~{}~{}~{}~{}~{}~{}~{}~{}~{}Q_{2}=(\bar{q}_{{}_{f}}\gamma^{5}q)(\bar{\chi}\chi),~{}~{}~{}~{}~{}~{}~{}~{}~{}Q_{3}=(\bar{q}_{{}_{f}}q)(\bar{\chi}\gamma^{5}\chi),$
(2) $\displaystyle
Q_{4}=(\bar{q}_{{}_{f}}\gamma^{5}q)(\bar{\chi}\gamma^{5}\chi),~{}~{}~{}Q_{5}=(\bar{q}_{{}_{f}}\gamma_{\mu}q)(\bar{\chi}\gamma^{\mu}\gamma^{5}\chi),~{}~{}~{}Q_{6}=(\bar{q}_{{}_{f}}\gamma_{\mu}\gamma^{5}q)(\bar{\chi}\gamma^{\mu}\gamma^{5}\chi),$
$\displaystyle
Q_{7}=(\bar{q}_{{}_{f}}\gamma_{\mu}q)(\bar{\chi}\gamma^{\mu}\chi),~{}~{}~{}Q_{8}=(\bar{q}_{{}_{f}}\gamma_{\mu}\gamma^{5}q)(\bar{\chi}\gamma^{\mu}\chi),~{}~{}~{}Q_{9}=(\bar{q}_{{}_{f}}\sigma_{\mu\nu}q)(\bar{\chi}\sigma^{\mu\nu}\chi).$

Two points should be stressed. First, in a specific model, a mediator maybe
assumed to connect the SM quarks and the dark sector. In that case, not all
the operators will contribute. Second, the coupling constants of some
operators could be suppressed severely from a more fundamental point of view.
For example, as Ref. Kamenik and Smith (2012) pointed out, if one starts from
an effective Lagrangian which respects the SM gauge symmetry, then the
dimension-seven operators $(H^{\dagger}\bar{q}_{{}_{f}}q)(\bar{\chi}\chi)$ and
$(H^{\dagger}\bar{q}_{{}_{f}}\sigma_{\mu\nu}q)(\bar{\chi}\sigma^{\mu\nu}\chi)$
should be included. After electroweak symmetry breaking, they are reduced to
$Q_{1}$ and $Q_{9}$, respectively, whose coefficients are suppressed by an
additional factor $\frac{v}{\Lambda}$ with $v$ being the vacuum expectation
value of Higgs field.

The upper limits of the coupling constants in the effective Lagrangian can be
achieved by comparing the difference between theoretical predictions and the
experimental data. As the corresponding detection of the $B_{c}$ meson is
still missing, we cannot use the experimental data of $B_{c}$ meson to set
constraints directly. Instead, the allowed region of the coupling constants
can be obtained by considering the $B$ meson decay processes. These channels
are $B^{-}\to K^{-}(K^{\ast-})+\not{E}$ and
$B^{-}\to\pi^{-}(\rho^{-})+\not{E}$, which have the same vertex as that of the
$B_{c}$ meson decays. The upper bounds of the $B$ meson decays involving
missing energy are listed in the first column of Table LABEL:tab1. These
results are dependent on experimental accuracy. With more precise experimental
conditions, these results might be further compressed in the future. Although
we have cited the most stringent results so far, all of them leave room for
contributions from new physics Grygier  _et al._ (2017). The second column is
the theoretical predictions, and the third one is the extracted upper limits
for the decays involving the assumed particles. One notices that they are of
the same order as that of the SM background.

Table 1: The branching ratios (in units of $10^{-6}$) of $B$ decays involving missing energy. Experimental bound Chen  _et al._ (2007); Grygier  _et al._ (2017); Lai  _et al._ (2017) | SM prediction Kamenik and Smith (2009); Jeon  _et al._ (2006); Altmannshofer  _et al._ (2009); Bartsch  _et al._ (2009) | Invisible particles bound
---|---|---
${\rm BR}(B^{\pm}\to K^{\pm}\not{E})<14$ | ${\rm BR}(B^{\pm}\to K^{\pm}\nu\bar{\nu})=5.1\pm 0.8$ | ${\rm BR}(B^{\pm}\to K^{\pm}\chi\chi)<9.7$
${\rm BR}(B^{\pm}\to\pi^{\pm}\not{E})<14$ | ${\rm BR}(B^{\pm}\to\pi^{\pm}\nu\bar{\nu})=9.7\pm 2.1$ | ${\rm BR}(B^{\pm}\to\pi^{\pm}\chi\chi)<6.4$
${\rm BR}(B^{\pm}\to K^{*\pm}\not{E})<61$ | ${\rm BR}(B^{\pm}\to K^{*\pm}\nu\bar{\nu})=8.4\pm 1.4$ | ${\rm BR}(B^{\pm}\to K^{*\pm}\chi\chi)<54$
${\rm BR}(B^{\pm}\to\rho^{\pm}\not{E})<30$ | ${\rm BR}(B^{\pm}\to\rho^{\pm}\nu\bar{\nu})=0.49^{+0.61}_{-0.38}$ | ${\rm BR}(B^{\pm}\to\rho^{\pm}\chi\chi)<30$

###  II.1 $\chi$ is a Majorana fermion

We first consider the situation that the invisible particle is a Majorana
fermion. In such a case, the vector and tenor currents give no contribution,
namely, $\bar{\chi}\gamma^{\mu}\chi=0$ and $\bar{\chi}\sigma^{\mu\nu}\chi=0$
(these are not true for the Dirac fermion). For the $0^{-}\to 0^{-}$
transitions, only three operators give nonzero contribution. The effective
Lagrangian reads

$\displaystyle\mathcal{L}_{1}=g_{m1}(\bar{q}_{{}_{f}}q)(\bar{\chi}\chi)+g_{m3}(\bar{q}_{{}_{f}}q)(\bar{\chi}\gamma^{5}\chi)+g_{m5}(\bar{q}_{{}_{f}}\gamma_{\mu}q)(\bar{\chi}\gamma^{\mu}\gamma^{5}\chi),$
(3)

where the subscript $m$ in $g_{mi}$ indicates that we are dealing with
Majorana fermions. The hadronic transition matrix elements can be expressed
as,

$\displaystyle\langle M_{f}^{-}|(\bar{q}_{{}_{f}}q)|M^{-}\rangle$
$\displaystyle=\frac{M^{2}-M_{f}^{2}}{m_{q}-m_{q_{{}_{f}}}}f_{0}(s),$ (4)
$\displaystyle\langle M_{f}^{-}|(\bar{q}_{{}_{f}}\gamma_{\mu}q)|M^{-}\rangle$
$\displaystyle=(P+P_{f})_{\mu}f_{+}(s)+(P-P_{f})_{\mu}\frac{M^{2}-M_{f}^{2}}{s}\big{[}f_{0}(s)-f_{+}(s)\big{]},$
$\displaystyle\langle
M_{f}^{-}|(\bar{q}_{{}_{f}}\sigma_{\mu\nu}q)|M^{-}\rangle$
$\displaystyle=i\big{[}P_{\mu}(P-P_{f})_{\nu}-P_{\nu}(P-P_{f})_{\mu}\big{]}\frac{2}{M+M_{f}}f_{T}(s),$

where $P$ and $P_{f}$ are the momenta of the initial or final mesons,
respectively; $m_{q}$ and $m_{q_{f}}$ are the masses of quarks; $s$ is defined
as $(P-P_{f})^{2}$; $f_{+}$, $f_{0}$, and $f_{T}$ are form factors. Here we
adopt the results of the LCSR method Ball and Zwicky (2005) to write the form
factors as,

$\displaystyle f_{0}(s)=\frac{r_{2}}{1-s/m_{fit}^{2}},$ (5) $\displaystyle
f_{+,T}^{K}(s)=\frac{r_{1}}{1-s/m_{R}^{2}}+\frac{r_{2}}{(1-s/m_{R}^{2})^{2}},$
$\displaystyle
f_{+,T}^{\pi}(s)=\frac{r_{1}}{1-s/m_{R}^{2}}+\frac{r_{2}}{1-s/m_{fit}^{2}},$

where the corresponding parameters $r_{1}$, $r_{2}$, $m_{R}$, and $m_{fit}$
are presented in Table LABEL:tab2.

Table 2: Parameters in the form factors of the $B\to\pi(K)$ processes Ball and Zwicky (2005). $F_{i}$ | $r_{1}$ | $r_{2}$ |  $m_{fit}^{2}$ (GeV2) |  $m_{R}$ (GeV)
---|---|---|---|---
$f^{K}_{0}$ | $0$ | $0.330$ | $37.46$ | $\cdots$
$f^{K}_{+}$ | $0.162$ | $0.173$ | $\cdots$ | $5.41$
$f^{K}_{T}$ | $0.161$ | $0.198$ | $\cdots$ | $5.41$
$f^{\pi}_{0}$ | $0$ | $0.258$ | $33.81$ | $\cdots$
$f^{\pi}_{+}$ | $0.744$ | $-0.486$ | $40.73$ | $5.32$
$f^{\pi}_{T}$ | $1.387$ | $-1.134$ | $32.22$ | $5.32$

By finishing the three-body phase space integral, we get the branching ratio

$\mathcal{BR}=\frac{1}{512\pi^{3}M^{3}\Omega\Gamma_{B^{-}}}\int\frac{ds}{s}\lambda^{1/2}(M^{2},s,M_{f}^{2})\lambda^{1/2}(s,m_{\chi}^{2},m_{\chi}^{2})\int
d\cos\theta\sum_{\lambda}|\mathcal{M}|^{2},$ (6)

where $\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2xy-2xz-2yz$ is the
K${\rm\ddot{a}}$llen function; $m_{\chi}$ is the mass of the invisible
particle; $\theta$ is the angle between the three-dimensional momenta
$\vec{P}_{\chi}$ and $\vec{P}_{f}$ in the center-of-momentum frame of the
invisible particles; $\Gamma_{B^{-}}$ is the total width of $B^{-}$ meson;
$\Omega=2$ originates from the final two invisible particles being identical
(Majarana fermion), and $\Omega=1$ when $\chi$ is the Dirac fermion. In the
square of the amplitude, there are interference terms which come from the
contribution of two different operators. These terms are proved to be zero
when the invisible particles are (pseudo)scalars Li  _et al._ (2019) .
However, they are not all zero when $\chi$ is a fermion, which makes the
calculations much more complicated. We will also calculate these terms, and
actually for some of them, the contribution cannot be ignored.

The partial width can be written as

$\displaystyle\Gamma=\int{dPS_{3}\big{(}\sum_{j}g_{mj}\mathcal{T}_{j}\big{)}^{\dagger}\big{(}\sum_{i}g_{mi}\mathcal{T}_{i}\big{)}}=\sum_{ij}g_{mj}g_{mi}\widetilde{\Gamma}_{ij},$
(7)

where we have taken $g_{ij}$ to be real for simplicity, and defined
$\widetilde{\Gamma}_{ij}=\int
dPS_{3}\mathcal{T}^{\dagger}_{j}\mathcal{T}_{i}$, which is independent of the
effective coupling constants. Some interference terms are zero by themselves
or cancel each other out. The nonzero terms are $\widetilde{\Gamma}_{11}$,
$\widetilde{\Gamma}_{33}$, $\widetilde{\Gamma}_{55}$, and
$\widetilde{\Gamma}_{35}$. In Fig. 2, we plot them as functions of $m_{\chi}$.
The solid and dashed lines represent the noninterference and interference
terms, respectively. One can see that the two different channels have similar
results, because the final mesons $K$ and $\pi$ have the same quantum number
and small masses compared with that of the $B$ meson. The noninterference
terms decrease when $m_{\chi}$ gets larger, because the phase space gets
smaller. Detailed calculation shows that $\widetilde{\Gamma}_{11}$ and
$\widetilde{\Gamma}_{33}$ are proportional to $(p_{1}\cdot
p_{2}-m_{\chi}^{2})$ and $(p_{1}\cdot p_{2}+m_{\chi}^{2})$, respectively,
where $p_{1}$ and $p_{2}$ are the momenta of two final invisible particles. So
$\widetilde{\Gamma}_{11}$ is smaller than $\widetilde{\Gamma}_{33}$ except
when $m_{\chi}=0$. $\widetilde{\Gamma}_{55}$ is less than
$\widetilde{\Gamma}_{33}$ as they are related to different effective
operators. The interference terms $\widetilde{\Gamma}_{35}$ and its complex
conjugate $\widetilde{\Gamma}_{53}$ are numerically equal. One can see that
they are zero when $m_{\chi}=0$ as they are proportional to $m_{\chi}^{2}$,
which is quite different with $\Gamma_{ii}$.

(a) $B^{-}\to K^{-}\bar{\chi}\chi$

(b) $B^{-}\to\pi^{-}\bar{\chi}\chi$

Figure 2: $\tilde{\Gamma}_{ij}$ for $B\rightarrow K(\pi)\bar{\chi}\chi$ with
$\chi$ being a Majorana fermion.

The effective Lagrangian for the $0^{-}\to 1^{-}$ process has the form

$\displaystyle\mathcal{L}_{2}=g_{m2}(\bar{q}_{{}_{f}}\gamma^{5}q)(\bar{\chi}\chi)+g_{m4}(\bar{q}_{{}_{f}}\gamma^{5}q)(\bar{\chi}\gamma^{5}\chi)+g_{m5}(\bar{q}_{{}_{f}}\gamma_{\mu}q)(\bar{\chi}\gamma^{\mu}\gamma^{5}\chi)+g_{m6}(\bar{q}_{{}_{f}}\gamma_{\mu}\gamma^{5}q)(\bar{\chi}\gamma^{\mu}\gamma^{5}\chi).$
(8)

The hadronic transition matrix elements are parametrized by form factors
$A_{0}$, $A_{1}$, $A_{2}$, $A_{3}$ $V$, $T_{1}$, $T_{2}$, and $T_{3}$ Bharucha
_et al._ (2016); Isgur and Wise (1990); Aliev  _et al._ (2011),

$\displaystyle\langle M_{f}^{*-}|(\bar{q}_{{}_{f}}\gamma^{5}q)|M^{-}\rangle$
$\displaystyle=-i\big{[}\epsilon\cdot(P-P_{f})\big{]}\frac{2M_{f}}{m_{q}+m_{q_{{}_{f}}}}A_{0}(s),$
(9) $\displaystyle\langle
M_{f}^{*-}|(\bar{q}_{{}_{f}}\gamma_{\mu}\gamma^{5}q)|M^{-}\rangle$
$\displaystyle=i\bigg{\\{}\epsilon_{\mu}(M+M_{f})A_{1}(s)-(P+P_{f})_{\mu}\frac{\epsilon\cdot(P-P_{f})}{M+M_{f}}A_{2}(s)$
$\displaystyle~{}~{}~{}-(P-P_{f})_{\mu}\big{[}\epsilon\cdot(P-P_{f})\big{]}\frac{2M_{f}}{s}\big{[}A_{3}(s)-A_{0}(s)\big{]}\bigg{\\}},$
$\displaystyle\langle M_{f}^{*-}|(\bar{q}_{f}\gamma_{\mu}q)|M^{-}\rangle$
$\displaystyle=\varepsilon_{\mu\nu\rho\sigma}\epsilon^{\nu}P^{\rho}(P-P_{f})^{\sigma}\frac{2}{M+M_{f}}V(s),$
$\displaystyle\langle
M_{f}^{*-}|(\bar{q}_{{}_{f}}\sigma_{\mu\nu}q)|M^{-}\rangle,$
$\displaystyle=i\bigg{\\{}\varepsilon_{\mu\nu\rho\sigma}\epsilon^{\rho}(P+P_{f})^{\sigma}T_{1}(s)-\varepsilon_{\mu\nu\rho\sigma}\epsilon^{\rho}(P-P_{f})^{\sigma}\frac{M^{2}-M_{f}^{2}}{s}$
$\displaystyle~{}~{}~{}\times\big{[}T_{1}(s)-T_{2}(s)\big{]}-\varepsilon_{\mu\nu\rho\sigma}(P+P_{f})^{\rho}(P+P_{f})^{\sigma}\big{[}\epsilon\cdot(P-P_{f})\big{]}$
$\displaystyle~{}~{}~{}\times\big{\\{}\frac{1}{M^{2}-M_{f}^{2}}T_{3}(s)-\frac{1}{s}\big{[}T_{1}(s)-T_{2}(s)\big{]}\big{\\}}\bigg{\\}},$

where $\epsilon$ is the polarization vector of the final meson, and the
$\varepsilon^{0123}=1$ convention is used.

The form factors are parametrized by Bharucha  _et al._ (2016),

$F_{i}(s)=P_{i}(s)\sum_{k}\alpha_{k}^{i}\,\left[z(s)-z(0)\right]^{k},$ (10)

where the pole structure is $P_{i}(s)=(1-s/m_{R,i}^{2})^{-1}$; $F_{1}$,
$F_{2}$, $F_{3}$, $F_{4}$, $F_{5}$, $F_{6}$, $F_{7}$, and $F_{8}$ represent
$A_{0}$, $A_{1}$, $A_{12}$, $A_{3}$, $V$, $T_{1}$, $T_{2}$, and $T_{23}$,
respectively. $A_{2}$ and $T_{3}$ can be deduced from the relations

$\displaystyle A_{3}(s)$
$\displaystyle=\frac{M+M_{f}}{2M_{f}}A_{1}(s)-\frac{M-M_{f}}{2M_{f}}A_{2}(s)$
(11) $\displaystyle A_{12}(s)$
$\displaystyle=\frac{(M+M_{f}){}^{2}(M^{2}-M_{f}^{2}-s)A_{1}(s)-\big{[}(M+M_{f})^{2}-s\big{]}\big{[}(M-M_{f})^{2}-s\big{]}A_{2}(s)}{16MM_{f}^{2}(M+M_{f})}$
$\displaystyle T_{23}(s)$
$\displaystyle=\frac{(M^{2}-M_{f}^{2})(M^{2}+3M_{f}^{2}-s)T_{2}(s)-\big{[}(M+M_{f})^{2}-s\big{]}\big{[}(M-M_{f})^{2}-s\big{]}T_{3}(s)}{8MM_{f}^{2}(M-M_{f})}.$

$z(s)$ is defined as

$z(s)=\frac{\sqrt{s_{+}-s}-\sqrt{s_{+}-s_{0}}}{\sqrt{s_{+}-s}+\sqrt{s_{+}-s_{0}}},$
(12)

where $s_{\pm}\equiv(M\pm M_{f})^{2}$ and $s_{0}\equiv
s_{+}(1-\sqrt{1-s_{-}/s_{+}})$. The related parameters are listed in Table
LABEL:tab3.

Table 3: Parameters in the form factors of the $B\to\rho(K^{*})$ processes with $k_{\rm max}=2$ Bharucha  _et al._ (2016). $F_{i}$ | $B\to K^{*}$ |  $m_{R,i}^{b\to s}/$GeV | $B\to\rho$ |  $m_{R,i}^{b\to d}/$GeV
---|---|---|---|---
$\alpha_{0}^{A_{0}}$ | $0.36\pm 0.05$ |  | $0.36\pm 0.04$ |
$\alpha_{1}^{A_{0}}$ | $-1.04\pm 0.27$ | $5.366$ | $-0.83\pm 0.20$ | $5.279$
$\alpha_{2}^{A_{0}}$ | $1.12\pm 1.35$ |  | $1.33\pm 1.05$ |
$\alpha_{0}^{A_{1}}$ | $0.27\pm 0.03$ |  | $0.26\pm 0.03$ |
$\alpha_{1}^{A_{1}}$ | $0.30\pm 0.19$ | $5.829$ | $0.39\pm 0.14$ | $5.724$
$\alpha_{2}^{A_{1}}$ | $-0.11\pm 0.48$ |  | $0.16\pm 0.41$ |
$\alpha_{0}^{A_{12}}$ | $0.26\pm 0.03$ |  | $0.30\pm 0.03$ |
$\alpha_{1}^{A_{12}}$ | $0.60\pm 0.20$ | $5.829$ | $0.76\pm 0.20$ | $5.724$
$\alpha_{2}^{A_{12}}$ | $0.12\pm 0.84$ |  | $0.46\pm 0.76$ |
$\alpha_{0}^{V}$ | $0.34\pm 0.04$ |  | $0.33\pm 0.03$ |
$\alpha_{1}^{V}$ | $-1.05\pm 0.24$ | $5.415$ | $-0.86\pm 0.18$ | $5.325$
$\alpha_{2}^{V}$ | $2.37\pm 1.39$ |  | $1.80\pm 0.97$ |
$\alpha_{0}^{T_{1}}$ | $0.28\pm 0.03$ |  | $0.27\pm 0.03$ |
$\alpha_{1}^{T_{1}}$ | $-0.89\pm 0.19$ | $5.415$ | $-0.74\pm 0.14$ | $5.325$
$\alpha_{2}^{T_{1}}$ | $1.95\pm 1.10$ |  | $1.45\pm 0.77$ |
$\alpha_{0}^{T_{2}}$ | $0.28\pm 0.03$ |  | $0.27\pm 0.03$ |
$\alpha_{1}^{T_{2}}$ | $0.40\pm 0.18$ | $5.829$ | $0.47\pm 0.13$ | $5.724$
$\alpha_{2}^{T_{2}}$ | $0.36\pm 0.51$ |  | $0.58\pm 0.46$ |
$\alpha_{0}^{T_{23}}$ | $0.67\pm 0.08$ |  | $0.75\pm 0.08$ |
$\alpha_{1}^{T_{23}}$ | $1.48\pm 0.49$ | $5.829$ | $1.90\pm 0.43$ | $5.724$
$\alpha_{2}^{T_{23}}$ | $1.92\pm 1.96$ |  | $2.93\pm 1.81$ |

For $B\rightarrow K^{\ast}(\rho)\bar{\chi}\chi$, we plot
$\widetilde{\Gamma}_{ij}$ as functions of $m_{\chi}$ in Fig. 3. One notices
that $\widetilde{\Gamma}_{66}$ is larger than the other terms. There is only
one interference term $\widetilde{\Gamma}_{46}$ which is nonzero. Its
contribution is negative. $\widetilde{\Gamma}_{22}$ and
$\widetilde{\Gamma}_{55}$ are quite close to each other. For $B\rightarrow
K^{\ast}\bar{\chi}\chi$, these two terms are almost coincident.

(a) $B^{-}\to K^{*-}\bar{\chi}\chi$

(b) $B^{-}\to\rho^{-}\bar{\chi}\chi$

Figure 3: $\tilde{\Gamma}_{ij}$ for $B\rightarrow
K^{\ast}(\rho)\bar{\chi}\chi$ with $\chi$ being a Majorana fermion.

###  II.2 $\chi$ is a Dirac fermion

For the Majorana fermion, there is neither vector nor tensor interaction,
while for the Dirac fermion, these two kinds of interactions also give
contribution. When the invisible particle is a Dirac fermion, the effective
Lagrangian has more operators. For the $0^{-}\to 0^{-}$ transition, it can be
written as

$\displaystyle\mathcal{L}_{3}=$ $\displaystyle
g_{d1}(\bar{q}_{{}_{f}}q)(\bar{\chi}\chi)+g_{d3}(\bar{q}_{{}_{f}}q)(\bar{\chi}\gamma^{5}\chi)+g_{d5}(\bar{q}_{{}_{f}}\gamma_{\mu}q)(\bar{\chi}\gamma^{\mu}\gamma^{5}\chi)+g_{d7}(\bar{q}_{{}_{f}}\gamma_{\mu}q)(\bar{\chi}\gamma^{\mu}\chi)$
(13)
$\displaystyle+g_{d9}(\bar{q}_{{}_{f}}\sigma_{\mu\nu}q)(\bar{\chi}\sigma^{\mu\nu}\chi),$

where $g_{di}$s are the phenomenological coupling constants between the
invisible Dirac fermions and quarks. $\widetilde{\Gamma}_{ij}$s are presented
in Fig. 4. One notices that they are about a half of that in the Majorana case
where $\chi$ and its antiparticle $\bar{\chi}$ are identical. In Fig. 4, one
also notices that there are additional terms $\widetilde{\Gamma}_{77}$,
$\widetilde{\Gamma}_{99}$ and $\widetilde{\Gamma}_{79}$, which represent
vector and tensor currents, since they are not zero when $\chi$ is a Dirac
fermion. Like above, $\widetilde{\Gamma}_{55}$ and $\widetilde{\Gamma}_{77}$
have same value when $m_{\chi}=0$. The $\widetilde{\Gamma}_{99}$ term
increases first, then decreases to zero when the phase space gets less. The
interference term $\widetilde{\Gamma}_{79}$ has the same trend as
$\widetilde{\Gamma}_{25}$, for it is proportional to $m_{\chi}^{2}$.

(a) $B^{-}\to K^{-}\bar{\chi}\chi$

(b) $B^{-}\to\pi^{-}\bar{\chi}\chi$

Figure 4: $\tilde{\Gamma}_{ij}$ for $B\rightarrow K(\pi)\bar{\chi}\chi$ with
$\chi$ being a Dirac fermion.

For $0^{-}\to 1^{-}$ processes , the effective Lagrangian can be written as

$\displaystyle\mathcal{L}_{4}=$ $\displaystyle
g_{d2}(\bar{q}_{{}_{f}}\gamma^{5}q)(\bar{\chi}\chi)+g_{d4}(\bar{q}_{{}_{f}}\gamma^{5}q)(\bar{\chi}\gamma^{5}\chi)+g_{d5}(\bar{q}_{{}_{f}}\gamma_{\mu}q)(\bar{\chi}\gamma^{\mu}\gamma^{5}\chi)+g_{d6}(\bar{q}_{{}_{f}}\gamma_{\mu}\gamma^{5}q)(\bar{\chi}\gamma^{\mu}\gamma^{5}\chi)$
(14)
$\displaystyle+g_{d7}(\bar{q}_{{}_{f}}\gamma_{\mu}q)(\bar{\chi}\gamma^{\mu}\chi)+g_{d8}(\bar{q}_{{}_{f}}\gamma_{\mu}\gamma^{5}q)(\bar{\chi}\gamma^{\mu}\chi)+g_{d9}(\bar{q}_{{}_{f}}\sigma_{\mu\nu}q)(\bar{\chi}\sigma^{\mu\nu}\chi).$

In Fig. 5, we plot the nonzero $\widetilde{\Gamma}_{ij}$s as functions of
$m_{\chi}$. One can see that the $\widetilde{\Gamma}_{ii}$ terms are not equal
to zero even when $m_{\chi}$ takes zero. As before, the interference terms
begin from zero and end up with zero when $m_{\chi}$ increases. Comparing with
Fig. 3, one can see that the additional term $\widetilde{\Gamma}_{99}$ related
to $Q_{9}$ is larger than others. However, this does not mean that the tensor
current makes much larger contribution to the partial width, for it also
depends on $g_{f9}$, which might be suppressed compared with other coupling
constants.

(a) $B^{-}\to K^{*-}\bar{\chi}\chi$

(b) $B^{-}\to\rho^{-}\bar{\chi}\chi$

Figure 5: $\tilde{\Gamma}_{ij}$ for $B\rightarrow
K^{\ast}(\rho)\bar{\chi}\chi$ with $\chi$ being a Dirac fermion.

##  III The decay modes of the $B_{c}$ meson

In Sec. II, the LCSR method is adopted to calculate the hadronic transition
amplitude in the FCNC processes of $B$ meson, where the final meson is light.
While for the $B_{c}$ meson decay modes, both initial and final mesons are
heavy. Under these circumstances, the BS method is a good choice to calculate
the hadronic transition amplitude. In this method, we can safely make an
instantaneous approximation when solving the Bethe-Salpeter equation fulfilled
by the wave functions of the heavy mesons. Details of how to solve the
instantaneous BS equation can be found in Kim and Wang (2004); Wang (2006).
The hadronic transition matrix element has the form

$\displaystyle\langle h^{-}|\bar{q}_{1}\Gamma^{\xi}b|B_{c}^{-}\rangle$
$\displaystyle=\int\frac{d^{3}q}{(2\pi)^{3}}{\rm
Tr}\left[\frac{\not{P}}{M}\overline{\varphi}_{P_{f}}^{++}(\vec{q}_{f})\Gamma^{\xi}\varphi_{P}^{++}(\vec{q})\right],$
(15)

where $\varphi^{++}_{P}$ and $\varphi^{++}_{P_{f}}$ are the wave functions of
the initial and final mesons, respectively; $\vec{q}$ and $\vec{q}_{{}_{f}}$
are the relative momentum of the quark and antiquark in the initial and final
mesons, respectively. In the Standard Model, $\nu$ and $\bar{\nu}$ lead to the
missing energy in the decay processes $B_{c}^{-}\to D_{s}^{(*)-}+\not{E}$ and
$B_{c}^{-}\to B^{(*)-}+\not{E}$. The branching ratio of former channels is of
the order of $10^{-7}\sim 10^{-6}$, while for the later ones, it is of the
order of $10^{-15}\sim 10^{-14}$. The exact results can be found in our
previous paper Li  _et al._ (2019).

###  III.1 $\chi$ is a Majorana fermion

The decay processes of $B_{c}$ meson to Majorana fermions are also described
by the effective Lagrangians in Eqs. (3) and (8). Using Eq. (15), we get
$\widetilde{\Gamma}_{ij}$s as functions of $m_{\chi}$, which are plotted in
Fig. 6. Although a different method is used to parametrize the form factors,
the results of $B_{c}\rightarrow D_{(s)}$ and $B_{c}\rightarrow
D^{\ast}_{(s)}$ are quite similar to those in Figs. 2 and 3, respectively,
because the main difference of these channels comes from the different
spectator quarks. We also consider the processes $B_{c}\rightarrow
B^{(\ast)}\bar{\chi}\chi$. One notices that $\widetilde{\Gamma}_{ij}$s of the
$c\to u$ processes are 2 orders of magnitude less than these of the $b\to
d(s)$ processes. This is because the phase space of the former channel is less
than that of the later one. The $\widetilde{\Gamma}_{55}$ term is smaller than
$\widetilde{\Gamma}_{33}$ in $b$ decays while it is lager in $c$ decays, which
means this operator is less sensitive to the phase space.

(a) $B_{c}^{-}\to D_{s}^{-}\bar{\chi}\chi$

(b) $B_{c}^{-}\to D_{s}^{*-}\bar{\chi}\chi$

(c) $B_{c}^{-}\to D^{-}\bar{\chi}\chi$

(d) $B_{c}^{-}\to D^{*-}\bar{\chi}\chi$

(e) $B_{c}^{-}\to B^{-}\bar{\chi}\chi$

(f) $B_{c}^{-}\to B^{*-}\bar{\chi}\chi$

Figure 6: $\tilde{\Gamma}_{ij}$ for $B_{c}\rightarrow
h^{(\ast)}\bar{\chi}\chi$ with $\chi$ being a Majorana fermion.

The next step is to set constraints for the coupling constants and calculate
the upper limits for the branching fractions of $B_{c}$ decays. In Sec.II, we
obtained $\widetilde{\Gamma}_{ij}$ for the decay processes of $B$ meson.
Considering the upper limits of the branching fractions of such channels, we
can extract the allowed parameter space for the effective coupling constants.
Here we use two different ways to make the calculation. First, we assume just
one effective coupling constant is nonzero, and its upper bound can be easily
achieved. Of course, different operators will give different results. Second,
we will scan the whole parameter space spanned by all the coupling constants
under all the constraints. The first method sets maximum allowed region of
each constant, namely, $g_{fi}$s. By second method we divide each regions into
400 bins. The program runs multiple bins and selects the maximum value of the
branching ratio of $B_{c}$ meson when the selection of constants does not make
the $\mathcal{B}R$ of $B$ meson beyond the experimental upper limit.

With the effective coupling constants achieved above, we calculate the upper
limits for of the branching fractions of $B_{c}$ decays. The results are shown
in Fig. 7, where the dashed lines represent those calculated in the first way
and the solid line corresponds to that of the second way. One can see that the
results of two different ways do not coincide in most $m_{\chi}$ regions. The
difference comes from the contribution of the interference terms. For the
$B_{c}^{-}\to P$ processes, the three cases $ij=11$, $ij=33$, and $Total$ have
the same value when $m_{\chi}=0$. At some points, $ij=55$ coincides with
$Total$. For the $B_{c}^{-}\to V\chi\chi$ processes, $ij=66$ coincides with
$Total$ when $m_{\chi}=0$. One notices that as $m_{\chi}$ increases, the
branching ratios (for $Total$) first increase slowly and then decrease
rapidly. This is a result of the competition between the phase space and the
effective coupling constants. The upper limits of the branching fractions of
$B_{c}\to P\chi\chi$ are 1 order of magnitude less than those of $B_{c}^{-}\to
V\chi\chi$, which is mainly due to different experimental bounds in Table
LABEL:tab1. The upper limits of branching ratios will be scaled down, as more
precise experimental limits are obtained in the future.

(a) $B_{c}^{-}\to D_{s}^{-}\bar{\chi}\chi$

(b) $B_{c}^{-}\to D^{-}\bar{\chi}\chi$

(c) $B_{c}^{-}\to D_{s}^{*-}\bar{\chi}\chi$

(d) $B_{c}^{-}\to D^{*-}\bar{\chi}\chi$

Figure 7: The upper limits of branching ratios of $B_{c}$ decays to Majorana
fermions.

In Fig. 8 we present the differential branching fractions as functions of $s$
which is defined as $s=(P-P_{f})^{2}$. As examples, three cases with
$m_{\chi}=0$ GeV, $0.25(M-M_{f})$, and $0.4(M-M_{f})$, respectively, are
considered. For comparison, the SM background with $\bar{\nu}\nu$ emission is
also plotted as blue dashed lines, which are less than those of the
$\bar{\chi}\chi$ emission channels in most regions of $s$. The left starting
point of the curves is the lower bound of $s$, which is determined by the mass
of $\chi$. The position of peaks of the distribution curves is almost
independent of $m_{\chi}$, which is at the region $s=16\sim 18$ $\rm GeV^{2}$.
We can see that the peak value gets larger as $m_{\chi}$ increases, because
the branching ratio increases with $m_{\chi}$ (until reaches its maximum value
around $m_{\chi}=1.5$ GeV) .

(a) $B_{c}^{-}\to D_{s}^{-}\bar{\chi}\chi$

(b) $B_{c}^{-}\to D^{-}\bar{\chi}\chi$

(c) $B_{c}^{-}\to D_{s}^{*-}\bar{\chi}\chi$

(d) $B_{c}^{-}\to D^{*-}\bar{\chi}\chi$

Figure 8: The differential branching ratios of $B_{c}$ decays to Majorana
fermions.

###  III.2 $\chi$ is a Dirac fermion

The similar analysis can also be applied to the Dirac fermions. The effective
Lagrangians take the same forms as those in Eqs. (13) and (14). In Fig. 9, we
plot $\widetilde{\Gamma}_{ij}$ as the function of $m_{\chi}$, which is about
half of the corresponding one in the Majorana case. One can see there are
several additional terms $\widetilde{\Gamma}_{77}$, $\widetilde{\Gamma}_{88}$,
and $\widetilde{\Gamma}_{99}$ which do not exist in the Majorana case. The
effective coupling constants obtained by comparing with the experimental
results are used to find the maximum values of the branching fractions which
are plotted in Fig. 10. In Figs. 10(a) and 10(b), we give the results of the
${B_{c}}\to P\bar{\chi}\chi$ processes. If we only consider the contribution
of $O_{7}$ or $O_{9}$, one can see the upper bound of the branching ratios,
which are labeled by $ij=77$ and $ij=99$, respectively, are less than those
resulted by other operators. This means that they do not affect the maximum
branch fractions obtained by the second way, namely considering the operators
altogether. This leads to the result that the upper limits of the branching
ratios of such channels are the same as those in the Majorana case.
Correspondingly, the differential branching fractions of two cases are also
the same with each other.

(a) $B_{c}^{-}\to D_{s}^{-}\bar{\chi}\chi$

(b) $B_{c}^{-}\to D_{s}^{*-}\bar{\chi}\chi$

(c) $B_{c}^{-}\to D^{-}\bar{\chi}\chi$

(d) $B_{c}^{-}\to D^{*-}\bar{\chi}\chi$

(e) $B_{c}^{-}\to B^{-}\bar{\chi}\chi$

(f) $B_{c}^{-}\to B^{*-}\bar{\chi}\chi$

Figure 9: $\tilde{\Gamma}_{ij}$ for $B_{c}\rightarrow
h^{(\ast)}\bar{\chi}\chi$ with $\chi$ being a Dirac fermion.

(a) $B_{c}^{-}\to D_{s}^{-}\bar{\chi}\chi$

(b) $B_{c}^{-}\to D^{-}\bar{\chi}\chi$

(c) $B_{c}^{-}\to D_{s}^{*-}\bar{\chi}\chi$

(d) $B_{c}^{-}\to D^{*-}\bar{\chi}\chi$

Figure 10: The upper limits of branching ratios of $B_{c}$ decays to Dirac
fermions.

In Figs. 10(c) and 10(d), the results of the ${B_{c}}\to V\bar{\chi}\chi$
processes are presented. The upper bound labeled by $Total$ (solid blue line)
is the same as that in Fig. 7(c) or 7(d) when $m_{\chi}$ is less than $1.2$ or
$1.3$ GeV. It becomes larger when $m_{\chi}$ continues to increase, as the
$\widetilde{\Gamma}_{77}$ term, which does not exist in the Majorana case,
will give the main contribution. So, in this range, there are some differences
between the upper bounds obtained in two cases (of course, if only $O_{7}$,
$O_{8}$, and $O_{9}$ give contribution, the Majorana case is not allowed). The
errors of $\widetilde{\Gamma}_{ij}$ are about $\pm 20\%$ from LCSR with
uncertainties of parameters in Table III. It will affect the upper limits of
the coupling constants. By varying the parameters in BS method by $\pm 5\%$,
the errors are from $\pm 7\%$ to $\pm 17\%$. Total errors of branching ratios
of $B_{c}$ meson decays in Figs. 7 and 10 are about $\pm 30\%$.

Correspondingly, the differential branching ratios, which are plotted in Fig.
11, should also show some differences with those in the Majorana case. When
$m_{\chi}=0.4(M-M_{f})$, the distribution curves have clearly different shapes
from those in Figs. 8(c) and 8(d). By varying the parameters in BS method by
$\pm 5\%$, the errors of distribution curves are less than $\pm 10\%$. Further
considering the uncertainties from LCSR in Table LABEL:tab3, the total errors
are about $\pm 30\%$. The distributions of Majorana and Dirac types still can
be distinguished in regions of $11-13.5~{}\rm{GeV}^{2}$ ($B_{c}\to D_{s}^{*}$)
or $12-14~{}\rm{GeV}^{2}$ ($B_{c}\to D^{*}$). This might provide a way to
distinguish between them. It is necessary to notice that these results are
upper limits when we assume that all operators contribute at the same time. If
only a few operators contribute to this process, for example, $Q_{7}$,
$Q_{8}$, and $Q_{9}$, the distinction between Majorana and Dirac will become
very obvious. In such cases, only the Dirac type final fermions are allowed.

(a) $B_{c}^{-}\to D_{s}^{*-}\bar{\chi}\chi$

(b) $B_{c}^{-}\to D^{*-}\bar{\chi}\chi$

Figure 11: The differential branching ratios of $B_{c}$ decays with
$m_{\chi}=0.4~{}(M-M_{f})$. The dashed lines represent the errors come from BS
method. The dotted lines represent the errors come from both BS and LCSR.

##  IV Conclusion

We studied the FCNC processes of the $B$ meson decaying to the invisible
spin-1/2 fermions. Both the Majorana and Dirac cases were considered. The
effective Lagrangians were introduced to describe the coupling between
invisible particles and quarks. By comparing the theoretical predictions of
BR$(B\rightarrow D_{(s)}^{(\ast)}\bar{\nu}\nu)$ and the experimental upper
bounds for BR$(B\rightarrow D_{(s)}^{(\ast)}\not{E})$, we derived the
constraints of the effective coupling constants. By scanning the allowed
parameter space, we derived the upper limits of the branching fractions for
the similar processes of $B_{c}$ meson. When the final meson was a
pseudoscalar, the upper limits of the branching fractions was of the order of
$10^{-6}$, and when the final meson was vector, it was of the order of
$10^{-5}$. These results were much larger than those of the SM background. The
differential branching fractions of Majorana and Dirac invisible particles
were of different shapes when $m_{\chi}$ was larger than $1.2$ or $1.3$ GeV.
This could provide a way to distinguish between the Majorana type particle
from the Dirac one.

##  V Acknowledgments

This work was supported in part by the National Natural Science Foundation of
China under Grant No. 11575048. We also thank the HPC Studio at Physics
Department of Harbin Institute of Technology for access to computing resources
through<EMAIL_ADDRESS>
## References

  * Bernstein  _et al._ (1985) J. Bernstein, L. S. Brown, and G. Feinberg, Phys. Rev. D32, 3261 (1985).
  * Srednicki  _et al._ (1988) M. Srednicki, R. Watkins, and K. A. Olive, Nucl. Phys. B310, 693 (1988).
  * Izaguirre  _et al._ (2015) E. Izaguirre, G. Krnjaic, P. Schuster, and N. Toro, Phys. Rev. Lett. 115, 251301 (2015).
  * Bertone  _et al._ (2005) G. Bertone, D. Hooper, and J. Silk, Phys. Rept. 405, 279 (2005).
  * Komatsu  _et al._ (2009) E. Komatsu _et al._ (WMAP), Astrophys. J. Suppl. 180, 330 (2009).
  * Lee and Weinberg (1977) B. W. Lee and S. Weinberg, Phys. Rev. Lett. 39, 165 (1977).
  * Pospelov  _et al._ (2008) M. Pospelov, A. Ritz, and M. B. Voloshin, Phys. Lett. B662, 53 (2008).
  * Hooper and Zurek (2008) D. Hooper and K. M. Zurek, Phys. Rev. D77, 087302 (2008).
  * Aprile  _et al._ (2017) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 119, 181301 (2017).
  * Gligorov  _et al._ (2018) V. V. Gligorov, S. Knapen, M. Papucci, and D. J. Robinson, Phys. Rev. D97, 015023 (2018).
  * Abdallah  _et al._ (2014) J. Abdallah _et al._ , (2014), arXiv:1409.2893 [hep-ph] .
  * Belyaev  _et al._ (2019) A. Belyaev, E. Bertuzzo, C. Caniu Barros, O. Eboli, G. Grilli Di Cortona, F. Iocco, and A. Pukhov, Phys. Rev. D99, 015006 (2019).
  * Boehm and Fayet (2004) C. Boehm and P. Fayet, Nucl. Phys. B683, 219 (2004).
  * Kusenko (2009) A. Kusenko, Phys. Rept. 481, 1 (2009).
  * Hambye (2009) T. Hambye, JHEP 01, 028 (2009).
  * Li  _et al._ (2019) G. Li, T. Wang, Y. Jiang, X.-Z. Tan, and G.-L. Wang, JHEP 03, 028 (2019).
  * Drees and Hajkarim (2018) M. Drees and F. Hajkarim, JHEP 12, 042 (2018).
  * Yue and Chu (2018) C.-X. Yue and J.-P. Chu, Phys. Rev. D 98, 055012 (2018).
  * Djouadi  _et al._ (2012) A. Djouadi, O. Lebedev, Y. Mambrini, and J. Quevillon, Phys. Lett. B709, 65 (2012).
  * Djouadi  _et al._ (2013) A. Djouadi, A. Falkowski, Y. Mambrini, and J. Quevillon, Eur. Phys. J. C73, 2455 (2013).
  * Arcadi  _et al._ (2015) G. Arcadi, Y. Mambrini, and F. Richard, JCAP 1503, 018 (2015).
  * Hisano  _et al._ (2018) J. Hisano, R. Nagai, and N. Nagata, JHEP 12, 059 (2018).
  * Barman  _et al._ (2019) B. Barman, D. Borah, P. Ghosh, and A. K. Saha, JHEP 10, 275 (2019).
  * Restrepo  _et al._ (2019) D. Restrepo, A. Rivera, and W. Tangarife, Phys. Rev. D100, 035029 (2019).
  * Abe and Sato (2019) T. Abe and R. Sato, Phys. Rev. D99, 035012 (2019).
  * Fiaschi  _et al._ (2019) J. Fiaschi, M. Klasen, and S. May, JHEP 05, 015 (2019).
  * Matsumoto  _et al._ (2019) S. Matsumoto, Y.-L. S. Tsai, and P.-Y. Tseng, JHEP 07, 050 (2019).
  * Yang (2016) K.-C. Yang, Phys. Rev. D94, 035028 (2016).
  * Chala  _et al._ (2015) M. Chala, F. Kahlhoefer, M. McCullough, G. Nardini, and K. Schmidt-Hoberg, JHEP 07, 089 (2015).
  * Bird  _et al._ (2004) C. Bird, P. Jackson, R. V. Kowalewski, and M. Pospelov, Phys. Rev. Lett. 93, 201803 (2004).
  * Bird  _et al._ (2006) C. Bird, R. V. Kowalewski, and M. Pospelov, Mod. Phys. Lett. A21, 457 (2006).
  * Badin and Petrov (2010) A. Badin and A. A. Petrov, Phys. Rev. D82, 034005 (2010).
  * Gninenko and Krasnikov (2015) S. N. Gninenko and N. V. Krasnikov, Phys. Rev. D92, 034009 (2015).
  * Barducci  _et al._ (2018) D. Barducci, M. Fabbrichesi, and E. Gabrielli, Phys. Rev. D98, 035049 (2018).
  * Kamenik and Smith (2012) J. F. Kamenik and C. Smith, JHEP 03, 090 (2012).
  * Bertuzzo  _et al._ (2017) E. Bertuzzo, C. J. Caniu Barros, and G. Grilli di Cortona, JHEP 09, 116 (2017).
  * Aaltonen  _et al._ (2016) T. A. Aaltonen _et al._ (CDF), Phys. Rev. D 93, 052001 (2016).
  * Burdin (2016) S. Burdin (ATLAS), AIP Conf. Proc. 1735, 030003 (2016).
  * Berezhnoy  _et al._ (2019) A. Berezhnoy, I. Belov, A. Likhoded, and A. Luhinsky, Mod. Phys. Lett. A 34, 1950331 (2019).
  * Aaij  _et al._ (2019) R. Aaij _et al._ (LHCb), Phys. Rev. D 100, 112006 (2019).
  * Grygier  _et al._ (2017) J. Grygier _et al._ (Belle), Phys. Rev. D 96, 091101 (2017).
  * Chen  _et al._ (2007) K. F. Chen _et al._ (Belle), Phys. Rev. Lett. 99, 221802 (2007).
  * Lai  _et al._ (2017) Y. T. Lai _et al._ (Belle), Phys. Rev. D95, 011102 (2017).
  * Kamenik and Smith (2009) J. F. Kamenik and C. Smith, Phys. Lett. B680, 471 (2009).
  * Jeon  _et al._ (2006) J. H. Jeon, C. S. Kim, J. Lee, and C. Yu, Phys. Lett. B636, 270 (2006).
  * Altmannshofer  _et al._ (2009) W. Altmannshofer, A. J. Buras, D. M. Straub, and M. Wick, JHEP 04, 022 (2009), arXiv:0902.0160 [hep-ph] .
  * Bartsch  _et al._ (2009) M. Bartsch, M. Beylich, G. Buchalla, and D. N. Gao, JHEP 11, 011 (2009).
  * Ball and Zwicky (2005) P. Ball and R. Zwicky, Phys. Rev. D71, 014015 (2005).
  * Bharucha  _et al._ (2016) A. Bharucha, D. M. Straub, and R. Zwicky, JHEP 08, 098 (2016).
  * Isgur and Wise (1990) N. Isgur and M. B. Wise, Phys. Rev. D42, 2388 (1990).
  * Aliev  _et al._ (2011) T. M. Aliev, M. Savci, and K.-C. Yang, Phys. Lett. B700, 55 (2011).
  * Kim and Wang (2004) C. Kim and G.-L. Wang, Phys. Lett. B 584, 285 (2004).
  * Wang (2006) G.-L. Wang, Phys. Lett. B 633, 492 (2006).


# Shape changes of a single hairy particle with mobile ligands at a liquid-
liquid interface

T. Staszewski  $|$ , M. Borówko  $|$

(Received July 5, 2023, in final form November 8, 2023)

###### Abstract

We investigate rearrangements of a single hairy particle at a liquid-liquid
interface using coarse-grained molecular dynamics simulations. We consider the
particles with the same (symmetrical interactions) and different (asymmetrical
interactions) affinity to the liquids. We show how ligand mobility affects the
behavior of the hairy particle at the liquid-liquid interface. We found that
such a hairy particle can take various shapes at the interface. For example, a
Janus-like snowman consisting of a segment cluster and a bare part of the
core, Saturn-like structures, and the core with a wide “plume” on one side. A
configuration of the particle at the interface is characterized by the
vertical displacement distance and the orientation of the particle relative to
the phase boundary. The selected descriptors are used to characterize the
shape of the segment cloud. We found that the shape of a particle and its
localization at the interface can be determined by tuning the interactions
with the liquids.

Key words: hairy particles, particle-laden layers, molecular dynamics

###### Abstract

Ми досліджуємо перегрупування однієї волосистої частинки на межі поділу
‘‘рідина-рідина’’ з використанням огрубленої молекулярної динаміки.
Розглядаються частинки з однаковою (симетричні взаємодії) та різною
(асиметричні взаємодії) спорідненістю до рідин. Показано, як рухливість
ліганду впливає на поведінку волосистої частинки на межі поділу ‘‘рідина-
рідина’’. З’ясовано, що така волосиста частинка може приймати різні форми на
межі розділу: наприклад, подібного на Януса сніговика, що складається зі
скупчення сегментів і оголеної частини ядра; структур, подібних до Сатурна;
ядра з широким ‘‘шлейфом’’ з одного боку. Конфігурація частинки на межі
розділу характеризується відстанню вертикального зміщення та орієнтацією
частинки відносно межі розділу фаз. Вибрані дескриптори використовуються для
характеристики форми сегментної хмари. Виявлено, що форму частинки та її
локалізацію на межі розділу можна визначити шляхом регулювання її взаємодії з
рідинами.

Ключов слова: волосисті частинки, збагачені частинками шари, молекулярна
динаміка

##  1 Introduction

Polymer-tethered (hairy) nanoparticles have received much attention in
nanomaterials science  [1, 2, 3]. They are composed of a core and a layer of
polymer chains grafted by one end on the core’s surface. Such hybrid particles
naturally combine the features of both ligands and cores. Their properties can
be tuned by altering the chemical nature, grafting density, chain length, and
grafting architecture of the polymer layer, or by changing the nanoparticle
size and shape. Due to their unique properties, hairy particles have a broad
range of applications, such as stabilization of Pickering emulsions,
production of composites, sensors, drug delivery systems, pollution removal,
etc.  [1, 2, 3].

Most of the research focused on the modelling of the morphology of polymer
layers. To solve this problem various theoretical approaches were used [4, 5,
6, 7, 8, 9]. Ohno et al. [4] proposed the mean-field theory of the polymer-
tethered spherical particles of different sizes which is an extension of the
Daud-Cotton method developed for star polymers [5]. The self-consistent field
model and the scaling theory were also used to study configurations of chains
tethered on spherical particles [6, 7]. Moreover, Lo Verso et al. [8] used the
density functional theory to study polymers end-grafted to spherical
nanoparticles under good solvent conditions. However, Ginzburg [9] used a
self-consistent field-density functional theory to investigate the self-
assembly of hairy particles. The internal structure of hairy particles was
also studied using molecular simulations [10, 11, 12, 13, 14, 15, 16, 17].
Dong and Zhou [10] investigated the responsive behavior of hairy particles in
different solvents. They considered the particles modified with block
copolymers or mixed polymer brushes and found various structures depending on
the type of ligands and solvents. The reconfiguration inside the polymer
canopy can result in the formation of patchy nanoparticles [14, 15]. Molecular
simulations also showed that the adsorption of small particles on ligands can
change the structure of the polymer shell [16, 17].  Moreover, the variability
of the internal structure of ligand-tethered particles considerably affects
their assembly. Numerous studies have focused on the self-assembly of hairy
particles in bulk systems, as these particles are promising building blocks
for the production of novel materials [18, 19, 20, 21, 22].

All theoretical studies have clearly shown that in bulk systems hairy
particles change their shapes in response to the environment. Their shape-
flexibility depends on the nature of ligands, solvents, grafting density, and
temperature. These results are a signpost for studying hairy particles at
fluid-solid and fluid-fluid interfaces.

The materials containing nanoparticles can release them into the environment.
Therefore, it is necessary to develop the methods of removing such pollutants.
The adsorption on solid surfaces was considered the most effective and safe
procedure for the removal of nanoparticles [3]. The amount of adsorption on
solids depends on the size and shape of the particles. Therefore, changes in
the shape of hairy particles on the solid surfaces and inside the pores were
investigated [23, 24, 25, 26, 27].

The behavior of the hairy particles at fluid-fluid interfaces was also
extensively studied both from experimental and theoretical perspectives [28,
29, 30, 31, 32, 33, 34]. Numerous medical applications of hairy particles
require a deep understanding of the mechanism of their attachment at
interfaces and the transfer to cells [2]. It has been well established that
surface-active moieties (e.g., surfactant, (bio)polymers, and micro- and
nanoparticles) tend to accumulate at fluid-fluid interfaces. The adsorption is
driven by minimization of the Gibbs free energy owing to a decrease in the
contact area between the two fluids. The energies associated with the
attachment of particles to fluid interfaces depend on the chemical nature,
size, roughness, and wettability of the particles characterized by their
contact angles [35, 28, 29]. Due to large sizes of nanoparticles, their
attatchement enegies at fluid-fluid interfaces are much higher than those of
standard amphiphiles. For this reason, the adsorption of the particles at the
droplet-continuous phase interface is essential for the stabilization of these
systems [3, 35, 36]. The properties of hairy particles can be tuned in a wide
range by changing the nature of polymer coating. Thus, they can be considered
as potentially highly efficient emulsifiers. The morphology of polymer
coatings plays a significant role in these processes. At liquid-liquid
interfaces, the chains on the two sides of the interface can adopt different
configurations depending on the nature of the ligands, the quality of both
solvents, the grafting density, and other parameters [29].

The structure of individual hairy particles adsorbed at the interface was
extensively investigated by molecular simulations [31, 32, 33, 34]. Tay and
Bresme simulated the alkylthiol passivated gold nanocrystals trapped at the
air-water interface and showed that the shape of the hairy particles was
strongly perturbed by the interface [30]. Quan et al. [32] studied the
internal structure of gold nanoparticles modified with amphiphilic brushes at
the oil-water interface. They discussed the influence of grafting architecture
(diblock, mixed, and Janus brush-grafted particles) and the hydrophilicity of
polymers on particle morphology. Quite recently, Tang et al. [33] studied the
hairy particles on the interface between two mutually immiscible polymers.
They concluded that the free energy of rearrangement of the ligands
significantly contributes to the total change in free energy upon adsorption
of the particle. The Brownian diffusion of various hairy particles at the
liquid-liquid interface was also simulated and a good agreement between
experimental and simulation results was obtained only when the flexibility of
the particle shape had been taken into account [34].

The analysis of the behavior of individual hairy particles at liquid-liquid
phase boundaries can be a starting point for the study of their self-assembly.
Hairy particles can form either disordered monolayers or different two-
dimensional structures of a high degree of ordering at these interfaces [3,
28, 29]. Menath et. al [37] proved that the changes in the shape of hairy
particles at the interface considerably influence their self-assembly.

As has been shown above, the capability of hairy particles to change their
shapes in various environments is important for numerous processes. Although
the influence of many factors on the morphology of the polymer coatings of
hairy particles has been assessed, the mobility of ligands has rarely been
taken into account. It is well known, however, that the properties of hairy
particles depend on the type of grafting bonds [29, 38]. For example, the
thiol-gold bond is mobile, allowing thiolated ligands to move along the
nanoparticle surface, while polyelectrolyte brushes grown from the surface of
silica particles are irreversibly attached [29]. Most theoretical models
assumed that the chains were permanently anchored at certain points on the
core surface [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. Only a few theoretical
works have dealt with mobile ligands [27, 18, 39]. In such approaches, the
binding segments of ligands can freely move on the core surface but they are
incapable of breaking away from it. Ligand mobility can significantly
influence the behavior of hairy particles. For example, under certain
conditions, the chains can form clumps on the core [14, 15]. This, in turn,
determines the self-assembly of hairy particles. Molecular simulations showed
that hairy particles with fixed and mobile ligands can assemble in completely
different structures [18].

In this work, we use molecular dynamics to study individual hairy particles at
liquid-liquid interfaces. Our study is motivated by the well-known phenomenon
of nanoparticles migrating to liquid-liquid interfaces to reduce their surface
energy. We take advantage of this effect to trap hairy particles at the
interface between two mutually immiscible liquids. The particles that interact
similarly with both liquids would be expected to localize at the interfacial
plane, where they can occlude the largest possible area of the interface and
thereby provide the largest free energy benefit. As a consequence, the
particle adsorption lowers the surface tension and stabilizes the interface
[3, 35, 36, 28, 29]. Nanoparticles are strong emulsifiers [36]. In the case of
particles that differently interact with the liquids forming the phase
boundary, these effects are stronger. As it has been mentioned, the ability of
nanoparticles to stabilize the emulsions depends on the area of the occupied
interface [3, 35, 36]. Therefore, we discuss here the shape transformations of
hairy particles at the liquid-liquid interface. We consider the particles with
mobile ligands. In the case of such particles, changes in shape, and thus in
the area of the occupied phase boundary, should be more visible. The effects
of ligand mobility on the behavior of hairy particles at liquid-liquid
interfaces have not been investigated so far.

The rest of the work is organized as follows. In section 2, we describe the
model used and the simulation method and define the calculated observables.
Our results are presented and discussed in section 3. Finally, in section 4 we
summarize and draw our conclusions.

##  2 Methods

###  2.1 Model

We consider a coarse-grained model of the system that contains two immiscible
fluids, $W$, $O$, and hairy particles. A single hairy particle is modelled as
a spherical core with attached linear chains. The polymer chains, are treated
as Kremer-Grest bead-chains [38], in which polymer segments, typically a few
monomers long, are represented by beads of size $\sigma_{S}$ and mass $m_{S}$.
We assume that all ligands have the same length (each chain consists of $M$
segments) and the number of ligands attached to each core is $f$. The cores
are treated as soft spheres of diameter $\sigma_{C}$ and mass $m_{C}$. To
reduce the number of parameters, we assume that all segments and molecules of
both fluids have the same diameters and masses.

Adjacent beads along the chain representing bonded segments are connected to
each other using a finitely extensible nonlinear elastic (FENE) segment-
segment potential

$U_{\text{FENE}}=-\frac{k}{2}R_{0}^{2}\ln\left[1-\left(\frac{r}{R_{0}}\right)^{2}\right],$
(2.1)

where $r$ is the separation distance between the segments, $k$ is the spring
constant and $R_{0}$ is the maximum possible length of the spring.  We
introduce a typical value of the constant $k=30$ [38]. For the chain backbone,
we assume that $R_{0}=1.5\sigma_{S}$.

The same potential is used to model the binding between the core and the first
segments of the ligands. In this case, however, $R_{0}=1.5\sigma_{CS}$, where
$\sigma_{CS}=0.5(\sigma_{C}+\sigma_{S})$. This potential allows the ligands to
slide along the surface of the core.

The interactions between nonbonded segments, liquid molecules, and cores are
treated using a cut-and-shifted Lennard-Jones potential [40]

$u^{(ij)}=\left\\{\begin{array}[]{ll}4\varepsilon_{ij}\left[(\sigma_{ij}/r)^{12}-(\sigma_{ij}/r)^{6}\
\right]+\Delta u^{(ij)}(r),&\ \ \ r<r_{\text{cut}}^{(ij)},\\\ 0,&\ \ \ {\rm
otherwise,}\end{array}\right.$ (2.2)

where

$\Delta u^{(ij)}(r)=-(r-r_{\text{cut}}^{(ij)})\partial
u^{(ij)}(r_{\text{cut}}^{(ij)})/\partial r,$ (2.3)

and $r^{(ij)}_{\text{cut}}$ is the cutoff distance,
$r^{(ij)}_{\text{cut}}=0.5(\sigma_{i}+\sigma_{j})$ $(i,j=C,S,O,W)$,
$\varepsilon_{ij}$ denotes the parameter that characterizes interaction
strengths between spherical species $i$ and $j$. The indices $C,S,W,O$
correspond to the cores, segments, liquid $W$, and liquid $O$, respectively.
We use the cutoff distance to switch on or switch off the attractive
interactions. For the attractive interactions,
$r^{(ij)}_{\text{cut}}=2.5\sigma_{(ij)}$ while, for repulsive interactions,
$r^{(ij)}_{\text{cut}}=\sigma_{(ij)}$.

Standard units are used herein. The diameter of segments is the distance unit,
$\sigma=\sigma_{S}$ while the segments-segment energy parameter,
$\varepsilon=\varepsilon_{SS}$ is the energy unit. The mass of a single
segment is the mass unit, $m_{S}=m$. The basic unit of time is
$\tau=\sigma(\varepsilon/m)^{1/2}$. The gravity effect is assumed to be
negligible.

We use standard reduced quantities, namely reduced distances $l^{*}=l/\sigma$,
reduced energies $E^{*}=E/\varepsilon$, and the reduced temperature is
$T^{*}=k_{\text{B}}T/\varepsilon$, where $k_{\text{B}}$ is the Boltzmann
constant, $T$ is an absolute temperature while $l$ and $E$ denote a distance
and an energy, respectively. The temperature is kept at $T^{*}=1$.

We also define the reduced density of the $i$th component as
$\rho^{*}_{k}=\rho_{k}\sigma_{k}^{3}$, where $\rho=N_{k}/V$ is its number
density, $N_{k}$ is the total number of “atoms” $k$ and V is the volume of the
system.

###  2.2 Simulation details

Molecular dynamics simulations were carried out using the LAMMPS package [41,
42]. The simulation box was a cuboid of reduced dimensions equal to
$L_{x}^{*}=L_{y}^{*}=L_{z}^{*}=L^{*}=89$ along the axes $x$, $y$, and $z$,
respectively. Standard periodic boundary conditions in the $x$ and $y$
directions were implemented. In the $z$-direction, two repulsive Lennard-Jones
(12-6) walls were used to close the system. The size $L_{z}^{*}$ was large
enough to obtain two uniform bulk phases: $W$-rich and $O$-rich ones. A
velocity-Verlet algorithm with a time step of $\Delta t=0.002\tau$ was used
for integrating the equations of motion, and a Nosé-Hoover thermostat with a
time constant $\tau=10$ was used for regulating the system temperature.

Each particle was constructed by placing the appropriate number of chains on
the core surface. First, a hairy particle was placed at the plane parallel to
the $x$ and $y$ axes at $z=0$. Then, equal numbers of molecules of the fluids
$W$ and $O$ were inserted into the simulation box. During simulations,
particles can freely rotate, translate and change their conformations.

The system was equilibrated for at least $10^{7}$ time steps until its total
energy reached a constant level, at which it fluctuated around the mean value.
The production runs were for at least $10^{7}$ time steps. To improve the
accuracy, all averages are obtained from several independent simulations.

We carried the simulations for the following parameters of the particles:
$\sigma_{C}^{*}=6$, $m_{C}^{*}=216$, $M=20$, and $f=20$. Similar parameters
were used by Tang et al. [33]. To reduce the number of parameters, we assumed
that the majority of interactions between the species are repulsive. The
considered interaction parameters are presented in the next section.

During simulations, we collected the density profiles of cores and segments,
energies of interactions between all species, and total energy, and calculated
various characteristics of the systems. Our goal was to estimate the position
of the adsorbed particle, its orientation, and the shape of the polymer
canopy.

###  2.3 Calculated observables

To characterize the position of particles at the interface, we estimated the
average distance of a core center from the phase boundary, i.e., the so-called
displacement distance $h^{*}_{c}$. However, the shape of a segment cloud
surrounding the core is characterized by the metrics generated by the gyration
tensor [43]

$G_{\alpha\beta}=\frac{1}{N^{\prime}}\Big{\langle}\sum_{i=1}^{N^{\prime}}{(r_{i,\alpha}-r_{0,\alpha})(r_{i,\beta}-r_{0,\beta})}\Big{\rangle},$
(2.4)

where $r_{i,\alpha}$ and $r_{0,\alpha}$ are the $\alpha$ component
($\alpha,\beta=x,y,z$ ) of the $i$-th segment and the center of mass of the
segment cloud, respectively.

After diagonalization of the gyration tensor we obtain its eigenvalues
$\lambda_{i}$ ($i=1,2,3$), ordered such that
$\lambda_{1}\geqslant\lambda_{2}\geqslant\lambda_{3}$. Then, we calculate
three invariants $I_{1}=\lambda_{1}+\lambda_{2}+\lambda_{3}$,
$I_{2}=\lambda_{1}\lambda_{2}+\lambda_{1}\lambda_{3}+\lambda_{2}\lambda_{3}$,
and $I_{3}=\lambda_{1}\lambda_{3}\lambda_{3}$. Finally, these values are used
to define the shape descriptors of the segment cloud [43].

The radius of gyration of the segment cloud is given by

$R_{g}^{2}=I_{1}=\frac{1}{N^{\prime}}\Big{\langle}\sum_{i=1}^{N^{\prime}}{\textbf{r}_{i0}^{2}}\Big{\rangle},$
(2.5)

where $\textbf{r}_{i0}=\textbf{r}_{i}-\textbf{r}_{0}$, and $\textbf{r}_{i}$
and $\textbf{r}_{0}$ are positions of the $i$-th segment and the center of
mass, respectively, and $N^{\prime}=fM$.

We resolve the vectors $\textbf{r}_{i0}$ into components parallel to the axes
$x$, $y$, $z$, and calculate the corresponding radii of gyration labeled
$R_{g\alpha}^{2}$ ($\alpha=x,y,z$).

Moreover, we define the relative shape anisotropy (asphericity) $\kappa$, the
prolateness $S^{*}$, and the acylindricity $C^{*}$ as

$\kappa=1-3\langle I_{2}/I^{2}_{1}\rangle,$ (2.6)
$S^{*}=\langle(3\lambda_{1}-I_{1})(3\lambda_{2}-I_{1})(3\lambda_{3}-I_{1})/I^{3}_{1}\rangle,$
(2.7) $C^{*}=\langle(\lambda_{2}-\lambda_{3})/I_{1}\rangle,$ (2.8)

where $\langle...\rangle$ denotes an average over all configurations.

The relative shape anisotropy takes on the values between 0 (perfectly
spherical objects) and 1 (rigid rods). For a regular planar array,
$\kappa=0.25$ [43]. The prolateness takes on the values between $-0.25$ and 2.
For prolate shapes $S^{*}>0$ while negative values of $S^{*}$ correspond to
oblate shapes. In a perfectly prolate segment cloud, the shorter axes of the
ellipsoid are the same ($\lambda_{1}>\lambda_{2}=\lambda_{3}$). On the
contrary, for perfectly oblate clouds, the longer axes are identical
($\lambda_{1}=\lambda_{2}>\lambda_{3}$). The acylindricity is $C^{*}\geqslant
0$ and for a perfect cylindrical symmetry $C^{*}=0$. It should be stressed,
however, that also a perfect sphere has a vanishing acylindricity.

##  3 Results and discussion

The behavior of hairy particles at the liquid-liquid interface depends on the
strengths of interactions between all single entities: cores, segments, and
liquid molecules. Various sets of interaction parameters were considered in
simulations of hairy particles [16, 17].

To model a phase boundary between immiscible liquids, we assume that
interactions between molecules of the same liquids are attractive, while
cross-interactions interactions are repulsive. We set
$\varepsilon^{*}_{WW}=\varepsilon^{*}_{OO}=\varepsilon^{*}_{WO}=1$ [33]. For
attractive interactions between the same molecules ($OO$ and $WW$)
$r^{(ij)}_{\text{cut}}=2.5\sigma_{(ij)}$, while repulsive interactions between
different molecules ($OW$) $r^{(ij)}_{\text{cut}}=\sigma_{(ij)}$. Under the
considered conditions, a very well-pronounced interface between these liquids
is formed. The density profiles of both liquids drop rapidly to zero in the
interphase region (see figure 1b).

We also introduce the following assumptions [26, 27, 44]: (i) interactions of
the cores with the remaining “atoms” are purely repulsive, (ii)  interactions
between the same “atoms” are attractive ($\varepsilon^{*}_{ii}=1$).

We consider the particles interacting (i) “symmetrically” and (ii)
“asymmetrically” with both liquids. In the first case, ligands attract the
molecules $O$ with the same strength as with the molecules $W$,
$\varepsilon^{*}_{SO}=\varepsilon^{*}_{SW}$. For the second type of particles
interactions with molecules $W$ are assumed to be repulsive whereas the
strength of attractive interactions with the molecules $O$ are attractive. We
carried out simulations for
$\varepsilon^{*}_{SO}=0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6$.

Notice that different sets of energy parameters $\varepsilon^{*}_{SO}$ and
$\varepsilon^{*}_{SW}$ can correspond to two types of real systems: (i)
different particles at the same interface (e.g., oil-water) or, conversely,
(ii) the same particle at the phase boundary between different liquids (more
and less polar organic solvents).

Our goal was to analyze how the mobility of the ligands and the strengths of
interactions between the segments and both liquids influence the behavior of
hairy particles at the interface, i. e., their shape position, and
orientation. Examples of the equilibrium configurations are depicted using
OVITO [45].

Figure 1: (Colour online) Results for a hairy particle with purely repulsive
ligands. (a) Exemplary configuration. The black sphere represents the core;
green spheres correspond to the segments bound to the core; red spheres
represent other segments. Yellow and blue dots correspond to molecules $O$ and
$W$, respectively. (b) Density profiles of the core (black line), the chain
segments (red line), the profiles of liquids $O$ (orange line) and $W$ (navy
line), $k=C,S,O,W$. The abscissa is scaled logarithmically. The green dashed
line represents the location of the interface.

The behavior of the considered particles in bulk liquids depends on the
strength of interactions between the ligands and liquid molecules [4, 6, 7, 8,
10, 15]. If these interactions are repulsive or weakly attractive, the mobile
ligands avoiding contact with the fluid form a cluster on the core’s surface.
In this way, the particle transforms into a kind of Janus-like particle
consisting of a naked part of the core and a cloud of segments. As these
interactions become stronger, the liquid molecules penetrate the polymer
layer, so that the chains gradually move away from each other and a typical
core-shell [4, 6, 7, 8, 10] structure is formed. Below we discuss the behavior
of such transformable particles at the liquid-liquid interface.

Figure 1a presents the exemplary configuration of a hairy particle with the
segments repelling both liquids. In this case, a Janus-like structure is
formed. Due to its shape, we refer to this configuration as  “snowman-shaped
particle” [46]. Such a particle adsorbs at the interface and “isolates” both
liquids. In other words, the particle covers a part of the interface which
leads to a decrease of a number of highly nonprofitable contacts between
molecules $O$ and $W$. The particle lies parallel to the interface. In figure
1b, the density profiles along the $z$-axis are presented. The core oscillated
around the position at the interface ($z^{*}=0$). The segment profile has one
relatively low and wide peak. The distribution of segments is symmetrical with
respect to the phase boundary.

It should be emphasized here that the attachment energy of Janus particles is
several times higher than that for corresponding homogeneous particles [3, 28,
29]. In this case, the particles exhibit amphiphilic properties. The
localization of a particle at the interface gives an additional energetic
profit associated with different interactions of particular parts of the
particle with both liquids. Janus particles have a greater capability to
stabilize the emulsions.

Figure 2: (Colour online) Equilibrium configurations for symmetrically
interacting hairy particles (side views), for different values of energy
parameters $\varepsilon^{*}_{SO}=\varepsilon^{*}_{WO}$: (a) 0.2, (b) 0.4, (c)
0.6, (d) 0.8, (e) 1.0, (f) 1.2, (g) 1.4, and (h) 1.6. The black sphere
represents the core; green spheres correspond to segments bound to the core;
red spheres represent other segments. Yellow and blue dots correspond to
molecules $O$ and $W$, respectively.

Now, we turn to the discussion of the behavior of particles with symmetrical
attractive interactions with the liquids. We have collected the configurations
of such particles in figure 2 (side view) and figure 3 (top view). As
previously, for weak attractive interactions, “snowman-shaped Janus-like
particles”  are observed (parts a-c). When the strength of these interactions
increases, the cloud of segments gradually “spills” at the liquid surface.
Liquid molecules penetrate the brush and the chains straighten up and they
move away from each other. In this way, a Saturn-like configuration is formed
(see parts d and e). Further increasing $\varepsilon^{*}_{SO}$ causes partial
pulling of the chains into the depths of both liquids. The energy gain
resulting from the interaction of the chains with the liquids outweighs the
gain associated with the reduction of the contact surface between the liquids
(parts g and h).

Figure 3: (Colour online) Equilibrium configurations for symmetrically
interacting hairy particles (views from above) for different values of energy
parameters $\varepsilon^{*}_{SO}=\varepsilon^{*}_{WO}$: (a) 0.2, (b) 0.4, (c)
0.6, (d) 0.8, (e) 1.0, (f) 1.2, (g) 1.4, and (h) 1.6. The black sphere
represents the core; green spheres correspond to the segments bound to the
core; red spheres represent other segments. Yellow and blue dots correspond to
molecules $O$ and $W$, respectively.

To quantify these observations, we calculated the corresponding density
profiles (see figure 4). In all cases, one narrow and symmetrical peak at
$z^{*}=0$ was obtained for the core. However, the segment density profiles
significantly vary as the parameter $\varepsilon^{*}_{SO}$ increases.
Obviously, the segment profiles are symmetrical and have maxima at $z^{*}=0$.
Initially, the segment distributions have typical Gaussian shapes and become
narrower and slightly higher (parts a-d). This reflects the spilling of the
polymers on the interface. For stronger interactions, the profiles become
wider and their shape considerably changes (parts e-h). The segment density
sharply decreases around $z^{*}=0$, then this decline is slower and finally
rapidly falls to zero.

Figure 4: (Colour online) Density profiles of the core ($k=C$, black lines)
and the chain segments ($k=S$, red lines) for symmetrically interacting hairy
particles and different values of energy parameters
$\varepsilon^{*}_{SO}=\varepsilon^{*}_{WO}$: (a) 0.2, (b) 0.4, (c) 0.6, (d)
0.8, (e) 1.0, (f) 1.2, (g) 1.4, and (h) 1.6. The green dashed lines represent
the location of the interface. The abscissas are scaled logarithmically.

A situation of particles with asymmetrical interactions at the interface is
much different from the previous case. Examples of their equilibrium
configurations are presented in figure 5. The ligands repel molecules $W$ and
attract molecules $O$. At the interface, the particles form the above-
mentioned Janus-like  snowman-shaped structures. For weak interactions with
the liquid $O$, the ligands form a compact “drop” of segments (parts a-d). As
the interactions become stronger, the chains gradually straighten in the
preferred liquid ($O$), avoiding the phase boundary (parts e-g). On the one
side of the core, an unfurled plume is formed. Moreover, with an increasing
value of $\varepsilon^{*}_{SO}$, the core shifts towards the liquid $O$, and
for sufficiently strong interactions, the particle leaves the interface (h).

Figure 5: (Colour online) Equilibrium configurations for asymmetrically
interacting hairy particles (side views) and different values of the energy
parameter $\varepsilon^{*}_{SO}$: (a) 0.2, (b) 0.4, (c) 0.6, (d) 0.8, (e) 1.0,
(f) 1.2, (g) 1.4, and (h) 1.6. Interactions between segments and molecules $W$
are repulsive. The black sphere represents the core; green spheres correspond
to segments bound to the core; red spheres represent other segments. Yellow
and blue dots correspond to molecules $O$ and $W$, respectively. The angle,
$\theta$ characterizing the particle orientation, and the displacement
distance, $h^{*}_{C}$ are shown in (c) and (h), respectively.

The analysis of the density profiles plotted in figure 6 confirms these
insights. For weak interactions (parts a-c), the core density has a maximum at
$z^{*}=0$. This means that the core is pinned to the interface. A further
increase in $\varepsilon^{*}_{SO}$ causes its partial entering into phase $O$
(parts d-g), and a maximum is shifted towards the liquid $O$ (unpinned
configurations). In the case of weak interactions, the segment density
profiles are symmetrical peaks located at $z^{*}>0$ (from the $O$-liquid
side). On the contrary, for stronger interactions, these peaks are
asymmetrical, extended towards the liquid $O$. However, when the particle goes
to the bulk phase, the segment distribution has a form of a symmetrical and
wide peak (part h).

Figure 6: (Colour online) Density profiles of the core ($k=C$, black lines)
and the chain segments ($k=S$, red lines) for asymmetrically interacting hairy
particles and different values of the energy parameter $\varepsilon^{*}_{SO}$:
(a) 0.2, (b) 0.4, (c) 0.6, (d) 0.8, (e) 1.0, (f) 1.2, (g) 1.4, and (h) 1.6.
Interactions between segments and molecules $W$ are repulsive. The abscissas
are scaled logarithmically. The green dashed lines represent the location of
the interface.

To characterize the orientation of particles, we use the angle between the
$xy$-plane and the line connecting the mass centers of the clouds of segments
and the center of the core (the angle $\theta$ is shown in figure 5c) [3, 28,
29]. Another parameter characterizing a particle configuration at the
interface is the average displacement depth, i.e., the distance of the center
of the core from the phase boundary[3, 28, 29] (see figure 5h). These
quantities are presented in figure 7 for particles interacting symmetrically
(black lines) and asymmetrically (red lines) as a function of
$\varepsilon^{*}_{SO}$. The results obtained for repulsive ligands are plotted
at $\varepsilon^{*}_{SO}=0$. We treat these interactions as infinitely weak
attractive ones.

Figure 7: (Colour online) (a) Orientation angle, $\theta$ and (b) the average
position of the core (displacement distance), $h^{*}_{C}$, for symmetrical
interactions particles (black circles) and asymmetrically interacting
particles with the segments repulsive for $W$ molecules (red circles) as
functions of the energy parameter $\varepsilon^{*}_{SO}$. The green dashed
lines represent the location of the interface.

In the case of symmetrical interactions (figure 7a), the particle lies
parallel to the phase boundary to occupy as much area on it as possible
(compare figure 2). However, for asymmetrical interactions with the liquid,
the particle can be differently oriented with respect to the phase boundary.
As the attractive interactions with one liquid become stronger, the particle
gradually changes its orientation, from tilted (figure 5 a–g) to normal
(figure 5h). For symmetrically interacting particles, the $h^{*}_{c}$
oscillates around zero (figure 7b). When the interactions are asymmetrical and
comparable or stronger than the interactions with the preferred liquid, the
displacement depth quickly increases. Such a behavior of hairy particles
results from a complex interplay between the reduction of the area of the
liquid-liquid interface caused by its adsorption and interactions between all
species. This is altered by the reconfiguration of a polymer canopy built of
mobile ligands. Qualitatively similar behavior was reported for Janus
ellipsoids and Janus dumbbells at liquid-liquid interfaces [3].

Figure 8: (Colour online) The density of segments in the surface layer
$z^{*}_{C}\in{<}-0.5,0.5{>}$ for symmetrically interacting particles (black
circles) and asymmetrically interacting particles with the segments repulsive
for $W$ molecules (red circles) as functions of the energy parameter
$\varepsilon^{*}_{SO}$.

We now turn our attention to the question of how the shape transformations at
the interface are affected by an increase in the strength of $SO$
interactions. Firstly, we analyze the average density of segments in the
surface layer $z^{*}_{C}\in{<}-0.5,0.5{>}$ (see figure 8). The segments
located here block the access to the interface. In the case of symmetrical
interactions, this density increases to a maximum at $\varepsilon^{*}_{SO}=1$
and then decreases. Indeed, for $\varepsilon^{*}_{SO}>1$, a part of the chains
starts to unfold in the direction of the liquid $O$. Thus, the effect of the
covering of the interface is the strongest for $\varepsilon^{*}_{SO}=1$.
However, for asymmetrically interacting particles, the surface density quickly
decreases, and for $\varepsilon^{*}_{SO}=1$ it reaches zero. With increasing
$\varepsilon^{*}_{SO}$, such particles change their orientations and
adsorption height ($h^{*}_{z}$). The segments “escape” from the surface.

Figure 9: (Colour online) Shape metrics for symmetrically interacting
particles (black circles) and asymmetrically interacting particles with the
segments repulsive for $W$ molecules (red circles) as functions of the energy
parameter $\varepsilon^{*}_{SO}$. (a) The average squared radius of gyration
in directions parallel to the liquid-liquid interface; (b) the $z$th component
of the squared radius of gyration; (c) The total radius of gyrations, and (d)
the relative shape anisotropy.

Further corroboration for the reconfigurable nature of the hairy particles is
offered by an analysis of their geometrical characteristics. For this purpose,
we employ the shape metrics computed starting from the gyration tensor. In
figure 9, we plot the total squared radius of gyration, its components and the
relative shape anisotropy as a function of $\varepsilon^{*}_{SO}$. In all
cases, the components in $x$ and $y$ directions are almost the same. To
improve the statistics, we calculated their averages, $R^{2}_{gxy}$. The
averaged squared radii of gyration in directions parallel to the interface are
presented in part a. The $z$-th components of $R^{2}_{g}$ are shown in part b,
while part c depicts the total squared radii of gyration. Figure 9d shows how
the shape anisotropy changes with an increase in the strength of $SO$
interactions.

We begin with the analysis of the curves obtained for symmetrically
interacting particles. As $\varepsilon^{*}_{SO}$ increases, the squared radius
of gyration $R^{2}_{gxy}$ rises and reaches the plateau at
$\varepsilon^{*}_{SO}=1$ (Saturn-like configurations). However, the
$z$-component initially slightly increases, then decreases to the minimum at
$\varepsilon^{*}_{SO}\approx 1$ and quickly increases. This well reflects the
effect of flattening the cloud of segments on the phase boundary. Obviously,
the later rise of the $R^{2}_{z}$ corresponds to the penetration of the chains
into the preferred liquid. The total squared radius of gyration varies in the
same way as the squared radius of gyration in the directions $x$ and $y$. The
asphericity of symmetrically interacting particles rapidly increases from 0.07
(repulsive interactions) to approximately 0.25 for
$0.6\leqslant\varepsilon^{*}_{SO}\leqslant 1.0$ and then quickly decreases.
This result confirms the formation of quasi-two-dimensional structures of
segments for $\varepsilon^{*}_{SO}=1$. With a further increase in
$\varepsilon^{*}_{SO}$, the polymer canopy becomes more and more spherical.

Other relations are found for asymmetrically interacting particles. As the
energy parameter increases to $\varepsilon^{*}_{SO}=1.4$, the squared radii of
gyration $R^{2}_{gxy}$, the $R^{2}_{gz}$ and $R^{2}_{g}$ monotonously
increase. Initially, the chains remain coiled but for stronger interactions
the polymer layer considerably swells. When the particle detaches from the
interface ($\varepsilon^{*}_{SO}=1.6$), the $R^{2}_{gxy}$ slightly decreases,
while $R^{2}_{gz}$ sharply rises. It is interesting that the relative shape
anisotropy is relatively low ($\kappa<0.1$) and remains almost unchanged for
all particles attached to the interface. For $\varepsilon^{*}_{SO}=$, however,
the $\kappa$ drops to a value close to the limit value, 0. Thus, the relative
shape anisotropy identified all the observed structures as almost spherical
ones.

Figure 10: (Colour online) The prolateness (a) and the acylindricity (b) for
symmetrically interacting particles (black circles) and asymmetrically
interacting particles with the segments repulsive for $W$ molecules (red
circles) as functions of the energy parameter $\varepsilon^{*}_{SO}$.

In figure 10 we report the prolateness $S^{*}$ and acylindricity $C^{*}$ as a
function of the energy parameter $\varepsilon^{*}_{SO}$. In the case of
symmetrical interactions, for weak interactions
($\varepsilon^{*}_{SO}\leqslant 0.6$), the cloud of segments becomes more and
more prolate while it is oblate at greater values of $\varepsilon^{*}_{SO}$. A
qualitative and rapid change is observed between values
$\varepsilon^{*}_{SO}=0$ and $\varepsilon^{*}_{SO}=0.8$. For
$\varepsilon^{*}_{SO}=1$, the cloud is almost perfectly oblate, and the
prolateness reaches the minimum value $S^{*}=-0.2$. These changes, well
reflect the transformations from  snowman-shaped particles to Saturn-like
structures and later “standing up” of chains from the surface, as shown in
figures 1, 2, and 3. This is also visible at the plot $C^{*}$ vs.
$\varepsilon^{*}_{SO}$. The acylindricity increases from 0.05 to a maximum
value $C^{*}=0.45$ at $\varepsilon^{*}_{SO}=1$. In a Saturn-like particle, the
polymer canopy has the shape of a hockey puck and deviates considerably from a
perfect cylinder of elongated shape. For stronger interactions, acylindricity
slowly decreases.

For asymmetrically interacting particles considered here, the prolateness is
always close to zero. However, for $\varepsilon^{*}_{SO}\leqslant 0.2$, the
$S^{*}$ is positive but for stronger interactions $S^{*}$ is negative. The
prolateness reaches a minimum value for $\varepsilon^{*}_{SO}=1$. However, the
acylindricity increases from 0.05 to 0.25 (for $\varepsilon^{*}_{SO}=1$) and
drops to 0.11. In this case, the relative shape anisotropy is also close to
zero. These two metrics classified our clouds as spherical.

Finally, it should be stressed that in the case of particles with permanently
anchored ligands, their shape changes are less spectacular [32, 37]. Such
hairy particles also retain a core-shell structure at the liquid-liquid
interface. Depending on the interactions with the liquids, their coronas take
on the shape of increasingly elongated ellipsoids [37]. However, the formation
of Janus-like structures has not been reported so far.

Our study shows how the hairy particles with mobile ligands change the shape
when they are trapped at a liquid-liquid interface.

##  4 Conclusions

We used molecular dynamics simulations to study individual spherical particles
modified with mobile ligands at the interface between partially miscible
Lennard-Jones liquids ($O/W$). To limit the number of parameters, we assumed
that almost all interactions are repulsive except the interactions between the
same molecules of liquids and the ligand-liquid interactions.

We considered two types of particles: particles symmetrically and
asymmetrically interacting with both liquids. In the first case, interactions
of chains with molecules of both liquids are the same. For the second group of
particles, interactions of ligands with the liquid $W$ are repulsive while
their interactions change from repulsive to strongly attractive. We found that
particles which in the same way interact with both liquids are always pinned
to the interface but a particle shape is strongly affected by the energy
parameter $\varepsilon^{*}_{SO}$.

Adsorption of symmetrically interacting particles results only from the
coverage of a part of the phase boundary, which leads to the reduction of
energetically unfavorable contacts between particles of different liquids.
Such particles are always pinned to the interface but a particle shape is
strongly affected by the energy parameter $\varepsilon^{*}_{SO}$. For weak
interactions, particles take on the form of  snowman-shaped Janus-like
particles consisting of a part of the bare core and a “drop” of segments and
orient parallel to the phase boundary. As interactions with liquid become
stronger, a Saturn-like structure is formed. A further increase in
$\varepsilon^{*}_{SO}$ causes the chains to be drawn deeper into the preferred
phase and particles take on the form of a core with a plume.

If interactions of particles with liquids are different, the surface activity
of particles can be additionally altered. In this case, we observe various
equilibrium configurations. The center of the core can be located either at
the interface (pinned configuration) or slightly shifted into a preferred
liquid (unpinned configuration). An increase of the energy parameter
$\varepsilon^{*}_{SO}$ causes the shape transformation of particles. For weak
interactions, Janus-like structures pinned to the interface are found.
Initially, the  particles lie at the interface but with increasing strength of
interactions the particles align themselves at an angle to the surface.  The
particle orientation gradually changes from parallel, through tilted to
perpendicular. This is similar to the behavior of Janus ellipsoids and Janus
dumbbells at the liquid-liquid interface [3]. In the case of stronger
interactions, the chains straighten in the preferred liquid and “escape” from
the interface. Currently, the core slightly shifts towards the interior of one
of the phases. On the one side of the core, an unfurled plume is visible.

We also discuss how the interactions with liquids affect the metric properties
characterizing the shape of segment clouds, such as the radius of gyration,
its components in the Cartesian coordinates, the relative shape anisotropy,
the prolateness, and acylindricity. Our analysis shows that the chosen
quantities well reflect the observed structures.

As shown in numerous studies [1, 3, 33, 37], particles with permanently
anchored ligands also change shapes at the liquid-liquid interface. The
particles with symmetrical interactions with both liquids deform into oblate
ellipsoids to maximize the occluded area at the interface. In the case of
different interactions with liquids, the polymer coronas become asymmetrical.
However, such particles do not form Janus-like and Saturn-like structures. We
conclude that ligand mobility increases the capability of hairy particles to
change the internal morphology at the interface. This strengthens their
surface activity, i.e., their tendency to stay at the interface. Due to this,
such particles can be effectively used to stabilize the emulsions.

We believe that our approach contains salient physical characteristics to
capture the mechanisms of rearrangement of hairy particles trapped at liquid-
liquid interfaces. The study proves that ligand mobility strongly affects the
behavior of an individual hairy particle at a liquid-liquid interface. Our
work may be a starting point for investigation of the impact of ligand
mobility on the adsorption of hairy particles at the interface, their self-
assembly and their capability to stabilize Pickering emulsions.

## References

  * [1] Moffitt M. G., J. Phys. Chem. Lett., 2013, 4, No. 21, 3654–3666, doi:10.1021/jz401814s.
  * [2] Pengo P., Şologan M., Pasquato L., Guida F., Pacor S., Tossi A., Stellacci F., Marson D., Boccardo S., Pricl S., Posocco P., Eur. Biophys. J., 2017, 46, No. 8, 749–771, doi:10.1007/s00249-017-1250-6.
  * [3] Borówko M., Staszewski T., Int. J. Mol. Sci., 2023, 24, No. 5, doi:10.3390/ijms24054564.
  * [4] Ohno K., Morinaga T., Takeno S., Tsujii Y., Fukuda T., Macromolecules, 2006, 39, No. 3, 1245–1249, doi:10.1021/ma0521708.
  * [5] Daoud M., Cotton J. P., J. Phys. France, 1982, 43, No. 3, 531–538, doi:10.1051/jphys:01982004303053100.
  * [6] Dan N., Tirrell M., Macromolecules, 1992, 25, No. 11, 2890–2895, doi:10.1021/ma00037a016.
  * [7] Wijmans C. M., Zhulina E. B., Macromolecules, 1993, 26, No. 26, 7214–7224, doi:10.1021/ma00078a016.
  * [8] Lo Verso F., Egorov S. A., Milchev A., Binder K., J. Chem. Phys., 2010, 133, No. 18, 184901, doi:10.1063/1.3494902.
  * [9] Ginzburg V. V., Macromolecules, 2017, 50, No. 23, 9445–9455, doi:10.1021/acs.macromol.7b01922.
  * [10] Dong J., Zhou J., Macromol. Theory Simul., 2013, 22, No. 3, 174–186, doi:10.1002/mats.201200078.
  * [11] Salerno K. M., Ismail A. E., Lane J. M. D., Grest G. S., J. Chem. Phys., 2014, 140, No. 19, 194904, doi:10.1063/1.4874638.
  * [12] Bolintineanu D. S., Lane J. M. D., Grest G. S., Langmuir, 2014, 30, No. 37, 11075–11085, doi:10.1021/la502795z.
  * [13] Giri A. K., Spohr E., J. Phys. Chem. C, 2018, 122, No. 46, 26739–26747, doi:10.1021/acs.jpcc.8b08590.
  * [14] Choueiri R. M., Galati E., Thérien-Aubin H., Klinkova A., Larin E. M., Querejeta-Fernández A., Han L., Xin H. L., Gang O., Zhulina E. B., Rubinstein M., Kumacheva E., Nature, 2016, 538, No. 7623, 79–83, doi:10.1038/nature19089.
  * [15] Staszewski T., J. Phys. Chem. C, 2020, 124, No. 49, 27118–27129, doi:10.1021/acs.jpcc.0c07775.
  * [16] Staszewski T., Borówko M., Phys. Chem. Chem. Phys., 2020, 22, 8757–8767, doi:10.1039/C9CP06854F.
  * [17] Borówko M., Staszewski T., Int. J. Mol. Sci., 2021, 22, No. 16, 8810, doi:10.3390/ijms22168810.
  * [18] Borówko M., Rżysko W., Sokołowski S., Staszewski T., Soft Matter, 2018, 14, 3115–3126, doi:10.1039/C8SM00213D.
  * [19] Rżysko W., Staszewski T., Borówko M., Colloids Surf., 2019, 570, 499–509, doi:https://doi.org/10.1016/j.colsurfa.2019.03.046.
  * [20] Ilnytskyi J. M., Slyusarchuk A., Saphiannikova M., Macromolecules, 2016, 49, No. 23, 9272–9282, doi:10.1021/acs.macromol.6b01871.
  * [21] Ilnytskyi J. M., Slyusarchuk A., Sokołowski S., Soft Matter, 2018, 14, 3799–3810, doi:10.1039/C8SM00356D.
  * [22] Slyusarchuk A., Yaremchuk D., Lintuvuori J., Wilson M. R., Grenzer M., Sokołowski S., Ilnytskyi J., Liq. Cryst., 2023, 50, No. 1, 74–97, doi:10.1080/02678292.2023.2169872.
  * [23] Ventura Rosales I. E., Rovigatti L., Bianchi E., Likos C. N., Locatelli E., Nanoscale, 2020, 12, 21188–21197, doi:10.1039/D0NR05058J.
  * [24] Che J., Park K., Grabowski C. A., Jawaid A., Kelley J., Koerner H., Vaia R. A., Macromolecules, 2016, 49, No. 5, 1834–1847, doi:10.1021/acs.macromol.5b02722.
  * [25] Staszewski T., Borówko M., Boguta P., J. Phys. Chem. B, 2022, 126, No. 6, 1341–1351, doi:10.1021/acs.jpcb.1c10418.
  * [26] Borówko M., Staszewski T., Int. J. Mol. Sci., 2022, 23, No. 14, 7919, doi:10.3390/ijms23147919.
  * [27] Staszewski T., Borówko M., Condens. Matter Phys., 2022, 25, No. 3, 33604, doi:10.5488/cmp.25.33604.
  * [28] Guzmán E., Martínez-Pedrero F., Calero C., Maestro A., Ortega F., Rubio R. G., Adv. Colloid Interface Sci., 2022, 302, 102620, doi:10.1016/j.cis.2022.102620.
  * [29] Garbin V., Crocker J. C., Stebe K. J., J. Colloid Interface Sci., 2012, 387, No. 1, 1–11, doi:10.1016/j.jcis.2012.07.047.
  * [30] Tay K. A., Bresme F., J. Am. Chem. Soc., 2006, 128, No. 43, 14166–14175, doi:10.1021/ja061901w.
  * [31] Lane J. M. D., Grest G. S., Phys. Rev. Lett., 2010, 104, 235501, doi:10.1103/PhysRevLett.104.235501.
  * [32] Quan X., Peng C., Dong J., Zhou J., Soft Matter, 2016, 12, 3352–3359, doi:10.1039/C5SM02721G.
  * [33] Tang T.-Y., Zhou Y., Arya G., ACS Nano, 2019, 13, No. 4, 4111–4123, doi:10.1021/acsnano.8b08733.
  * [34] Wang D., Zhu Y.-L., Zhao Y., Li C. Y., Mukhopadhyay A., Sun Z.-Y., Koynov K., Butt H.-J., ACS Nano, 2020, 14, No. 8, 10095–10103, doi:10.1021/acsnano.0c03291.
  * [35] Maestro A., Guzmán E., Ortega F., Rubio R. G., Curr. Opin. Colloid Interface Sci., 2014, 19, No. 4, 355–367, doi:10.1016/j.cocis.2014.04.008.
  * [36] Low L. E., Siva S. P., Ho Y. K., Chan E. S., Tey B. T., Adv. Colloid Interface Sci., 2020, 277, 102117, doi:10.1016/j.cis.2020.102117.
  * [37] Menath J., Eatson J., Brilmayer R., Andrieu-Brunsen A., Buzza D. M. A., Vogel N., PNAS, 2021, 118, No. 52, e2113394118, doi:10.1073/pnas.2113394118.
  * [38] Kremer K., Grest G. S., J. Chem. Phys., 1990, 92, No. 8, 5057–5086, doi:10.1063/1.458541.
  * [39] Borówko M., Sokołowski S., Staszewski T., Pizio O., J. Chem. Phys., 2018, 148, No. 4, 044705, doi:10.1063/1.5010687.
  * [40] Toxvaerd S., Dyre J. C., J. Chem. Phys., 2011, 134, No. 8, 081102, doi:10.1063/1.3558787.
  * [41] Thompson A. P., Aktulga H. M., Berger R., Bolintineanu D. S., Brown W. M., Crozier P. S., in ’t Veld P. J., Kohlmeyer A., Moore S. G., Nguyen T. D., Shan R., Stevens M. J., Tranchida J., Trott C., Plimpton S. J., Comput. Phys. Commun., 2022, 271, 108171, doi:10.1016/j.cpc.2021.108171.
  * [42] Plimpton S., J. Comput. Phys., 1995, 117, No. 1, 1–19, doi:10.1006/jcph.1995.1039.
  * [43] Theodorou D. N., Suter U. W., Macromolecules, 1985, 18, No. 6, 1206–1214, doi:10.1021/ma00148a028.
  * [44] Borówko M., Staszewski T., Tomasik J., J. Phys. Chem. B, 2023, 127, No. 22, 5150–5161, doi:10.1021/acs.jpcb.3c01943.
  * [45] Stukowski A., Modell. Simul. Mater. Sci. Eng., 2009, 18, No. 1, 015012, doi:10.1088/0965-0393/18/1/015012.
  * [46] Kang C., Honciuc A., ACS Nano, 2018, 12, No. 4, 3741–3750, doi:10.1021/acsnano.8b00960.

Ukrainian 0

Зміни форми волосистої частинки з рухливими лігандами на межі розділу
‘‘рідина-рідина’’ [Т. Сташевскі, М. Борувко]Т. Сташевскі, М. Борувко

Факультет теоретичної фімії, Інститут хімічних наук, Університет ім. Марії
Склодовської-Кюрі, Люблін, Польща