# On Coupling Constraints in Linear Bilevel Optimization

Dorothee Henke, Henri Lefebvre, Martin Schmidt, Johannes Thürauf

(Date: August 28, 2024)

###### Abstract.

It is well-known that coupling constraints in linear bilevel optimization can
lead to disconnected feasible sets, which is not possible without coupling
constraints. However, there is no difference between linear bilevel problems
with and without coupling constraints w.r.t. their complexity-theoretical
hardness. In this note, we prove that, although there is a clear difference
between these two classes of problems in terms of their feasible sets, the
classes are equivalent on the level of optimal solutions. To this end, given a
general linear bilevel problem with coupling constraints, we derive a
respective problem without coupling constraints and prove that it has the same
optimal solutions (when projected back to the original variable space).

###### Key words and phrases:

Bilevel optimization, Coupling constraints

###### 2020 Mathematics Subject Classification:

90Cxx

##  1\. Introduction

The research interest in bilevel optimization problems increased significantly
over the last years and decades; see, e.g., [6] for a recent overview.
However, and although serious advances have been made both w.r.t. theoretical
aspects and algorithmic developments, there are still open questions. This is
even the case for linear bilevel optimization problems that we consider in
this note and that are given by111For the ease of presentation, we omit
stating dimensions of matrices and vectors.

$\displaystyle\min_{x\in X,y}\quad$ $\displaystyle c^{\top}x+d^{\top}y$ (1a)
s.t. $\displaystyle Ax+By\geq a,$ (1b) $\displaystyle y\in S(x),$ (1c)

where $S(x)$ is the set of optimal solutions to the $x$-parameterized lower-
level problem

$\min_{y}\quad f^{\top}y\quad\text{s.t.}\quad Cx+Dy\geq b.$ (2)

In particular, we consider the optimistic linear bilevel problem, i.e., the
leader is able to choose the $y$ that is the best w.r.t. the upper-level
objective function if there are multiple optimal solutions to the follower’s
problem. Moreover, Problem (1) contains coupling constraints in (1b), i.e.,
upper-level constraints that explicitly depend on the lower-level variables.
These coupling constraints are the main topic of this note. Instead, a bilevel
problem with $B=0$ does not have any coupling constraints. Throughout the
remainder of the paper, we assume that the bilevel problem (1) is solvable and
that $X$ is a given polyhedron.

For motivating our main research question, we briefly discuss coupling
constraints in the following both w.r.t. their impact on the geometry of the
bilevel feasible set and their impact on complexity.

###  1.1. Geometry of the Feasible Set

As an example, we consider the following linear bilevel problem taken from
[12]:

$\displaystyle\min_{x,y}\quad$ $\displaystyle F(x,y)=x+6y$ (3a) s.t.
$\displaystyle-x+5y\leq 12.5,$ (3b) $\displaystyle y\in S(x),$ (3c)

with $S(x)$ being the set of optimal solutions to the lower-level problem

$\displaystyle\min_{y}\quad$ $\displaystyle f(x,y)=-y$ (4a) s.t.
$\displaystyle 2x-y\geq 0,$ (4b) $\displaystyle-x-y\geq-6,$ (4c)
$\displaystyle-x+6y\geq-3,$ (4d) $\displaystyle x+3y\geq 3.$ (4e)

Both levels are linear optimization problems and all variables are continuous.

$x$$y$$1$$2$$3$$4$$5$$6$1234$f$$F$$(\frac{3}{7},\frac{6}{7})$ Figure 1.
Visualization of Problem (3); mainly taken from [12].

The problem is visualized in Figure 1. The gray area is the set of points that
satisfy the lower-level constraints. The points above the dashed line are
infeasible w.r.t. the coupling constraint (3b). Due to the optimization
direction of the follower, the green faces denote the bilevel feasible set.
Note that the two red faces are not part of the bilevel feasible set since,
for the respective $x$-values, the optimal replies $y$ by the follower violate
the coupling constraint.

Hence, the example shows that the usage of coupling constraints makes it
possible to model bilevel feasible sets that are disconnected. Theorem 3.3 in
[3] states that the feasible set of a linear bilevel problem with $B=0$ is
always connected, which means that disconnected sets can only be modeled using
coupling constraints—an aspect that gained some prominence since it allows to
model mixed-binary linear problems using purely continuous linear bilevel
models; see, e.g., Section 3 in [15] and Section 3.1 in [1]. Consequently, it
seems to be the case that having coupling constraints introduces larger
modeling capabilities.

###  1.2. Complexity

Already [11] showed that linear bilevel problems are NP-hard. An alternative
proof using a reduction from KNAPSACK can be found in [2]. Moreover, [10]
proved that linear bilevel problems are strongly NP-hard by a reduction from
KERNEL; see Problem GT57 in [8]. In particular, they even showed that the
special case of min-max problems without coupling constraints is strongly NP-
hard. The same holds true for the reduction from 3-SAT shown in [13]. In [14]
it is further shown that even checking local optimality of a given point is
strongly NP-hard (again via 3-SAT). As before, the authors do not require
coupling constraints to achieve this hardness result; see the linear bilevel
problems in the proofs of Theorems 5.1 and 5.2 in [14]. Hence, coupling
constraints are not required to make linear bilevel problems strongly NP-hard.
Moreover, [4] recently showed that the decision versions of linear bilevel
problems with coupling constraints are in NP, which again implies that linear
bilevel problems with coupling constraints are not harder than those without.

###  1.3. Research Question

Combining the two last discussions, we can summarize the following two
conclusions:

  1. (i)

Using coupling constraints in linear bilevel problems allows for modeling a
richer class of feasible sets.

  2. (ii)

These stronger modeling capabilities do not result in any change of the
hardness of the resulting problem in terms of complexity theory.

Due to (ii), the question arises if we “really” increase modeling capabilities
by using coupling constraints. While this is the case w.r.t.
(dis-)connectedness of feasible sets, we prove that there is no difference on
the level of optimal solutions. To do this, in the next section, given a
linear bilevel problem with coupling constraints, we derive a linear bilevel
problem without coupling constraints that has the same set of optimal
solutions.

##  2\. Exact Penalization of Coupling Constraints

In this section, we show that the bilevel problem (1) with coupling
constraints can be reformulated as a bilevel problem without coupling
constraints. To do this, we show that the violation of the coupling
constraints (1b) can be exactly penalized in the objective function of the
leader. This result is surprising for at least two reasons. First, we just saw
that the feasible region of (1) may be disconnected and nonconvex. In this
case, standard Lagrangian duality theory is usually limited and not as strong
as in the convex case. Second, the resulting problem, i.e., after
penalization, has a smooth objective function and no coupling constraint.
Thus, it differs from classic exact penalty methods, which often require a
nonsmooth penalty function.

Our key idea is to reformulate Problem (1) in a way so that the follower
measures the violation of the coupling constraints directly in the lower-level
problem, while the leader enforces that this violation is zero. Doing so leads
to the following bilevel problem which contains a scalar and very simple
coupling constraint:

$\displaystyle\min_{x,y,\varepsilon}\quad$ $\displaystyle c^{\top}x+d^{\top}y$
(5a) s.t. $\displaystyle x\in X,$ (5b) $\displaystyle\varepsilon=0,$ (5c)
$\displaystyle(y,\varepsilon)\in\tilde{S}(x).$ (5d)

Here, $\tilde{S}(x)$ is the set of optimal solutions to the $x$-parameterized
lower-level problem

$\displaystyle\min_{y,\varepsilon}\quad$ $\displaystyle f^{\top}y$ (6a) s.t.
$\displaystyle Ax+By+\varepsilon e\geq a,$ (6b) $\displaystyle Cx+Dy\geq b,$
(6c) $\displaystyle\varepsilon\geq 0,$ (6d)

where $e$ is the vector of ones. Essentially, $\varepsilon$ is an additional
variable of the follower that measures the violation of the coupling
constraints. The newly introduced coupling constraint (5c) enforces that it
equals zero. Most importantly, we have $S(x)=\text{proj}_{y}(\tilde{S}(x))$
for all leader’s decisions $x$. In the next lemma, we show that Problem (5)
is, indeed, equivalent to Problem (1).

######  Lemma 2.1.

For every bilevel feasible point $(x,y)$ of Problem (1), the point $(x,y,0)$
is bilevel feasible for Problem (5) with the same objective value. For every
bilevel feasible point $(x,y,\varepsilon)$ of Problem (5), the point $(x,y)$
is bilevel feasible for Problem (1) with the same objective value.

###### Proof.

Let $(x,y)$ be a bilevel feasible point of Problem (1). Then, $(x,y,0)$
satisfies the upper- and lower-level constraints of Problem (5). The point
$(y,0)$ is optimal for the lower-level problem (6) since $\varepsilon$ is not
part of the lower-level objective function and $f^{\top}y$ is minimal
regarding Constraints (6c) due to the bilevel feasibility of $(x,y)$ for
Problem (1).

Let now $(x,y,\varepsilon)$ be a bilevel feasible point of Problem (5). Then,
the point $(x,y)$ satisfies the upper- and lower-level constraints of Problem
(1) because $\varepsilon=0$ holds. The point $y$ is optimal for the lower-
level problem (2) since every optimal solution $\tilde{y}$ to the
$x$-parameterized lower-level problem (2) can be extended to the feasible
point $(\tilde{y},\max\\{\left\lVert a-Ax-By\right\rVert_{\infty},0\\})$ of
the $x$-parameterized lower-level problem (6) with the same objective value.

We finally note that Problem (1) and Problem (5) have the same upper-level
objective functions, which proves the claim. ∎

We now penalize the single coupling constraint (5c) of Problem (5) to obtain a
bilevel problem without coupling constraints. Moreover, we show that there is
a polynomial-sized (in the bit-encoding length of the original problem’s data)
penalty parameter so that this formulation is equivalent in terms of optimal
solutions.

######  Theorem 2.2.

There is a polynomial-sized parameter $\kappa>0$ so that the bilevel problem
(without coupling constraints)

$\displaystyle\min_{x,y,\varepsilon}\quad$ $\displaystyle
c^{\top}x+d^{\top}y+\kappa\varepsilon$ (7a) s.t. $\displaystyle x\in X,\
(y,\varepsilon)\in\tilde{S}(x),$ (7b)

has the same set of optimal solutions as Problem (5). Again, $\tilde{S}(x)$ is
the set of optimal solutions to the $x$-parameterized lower-level problem (6).

The idea of the proof is as follows. First, we derive a single-level
reformulation of the bilevel problem (5), using the KKT conditions of the
follower’s problem (6). Second, we apply results from augmented Lagrangian
duality theory for mixed-integer linear problems to show that a polynomial-
sized exact penalization parameter exists. Finally, we show that the resulting
mixed-integer linear program is nothing but the KKT reformulation of Problem
(7).

###### Proof.

Since the lower-level problem (6) of Problem (5) is a linear program, we can
replace it with its KKT conditions [5], leading to

$\displaystyle\min_{x,y,\varepsilon}\quad$ $\displaystyle c^{\top}x+d^{\top}y$
(8a) s.t. $\displaystyle x\in X,\ \varepsilon=0,$ (8b) $\displaystyle
Ax+By+\varepsilon e\geq a,\ Cx+Dy\geq b,\ \varepsilon\geq 0,$ (8c)
$\displaystyle B^{\top}\lambda+D^{\top}\mu=f,\ e^{\top}\lambda+\eta=0,$ (8d)
$\displaystyle\lambda,\mu,\eta\geq 0,$ (8e)
$\displaystyle\lambda^{\top}(Ax+By+\varepsilon e-a)=0,\
\mu^{\top}(Cx+Dy-b)=0,\ \eta\varepsilon=0.$ (8f)

Using additional binary variables $z^{\lambda},z^{\mu},z^{\eta}$, and a
sufficiently large big-$M$ value, we can reformulate the complementarity
constraints (8f) as the mixed-integer linear constraints

$\displaystyle\lambda\leq(e-z^{\lambda})M,\ \mu\leq(e-z^{\mu})M,\
\eta\leq(1-z^{\eta})M,$ (9a) $\displaystyle Ax+By+\varepsilon e-a\leq
z^{\lambda}M,\ Cx+Dy-b\leq z^{\mu}M,\ \varepsilon\leq z^{\eta}M.$ (9b)

It is shown in [4] that a valid and polynomial-sized value for $M$ can be
computed in polynomial time. Thus, the resulting problem is a mixed-integer
linear program whose input data is polynomial-sized in the bit-encoding length
of the input data of Problem (5). We now penalize the constraint
$\varepsilon=0$ in the $\ell_{\infty}$-sense and obtain

$\displaystyle\min_{x,y,\varepsilon}\quad$ $\displaystyle
c^{\top}x+d^{\top}y+\kappa\varepsilon$ (10a) s.t. $\displaystyle x\in X,$
(10b) $\displaystyle Ax+By+\varepsilon e\geq a,\ Cx+Dy\geq b,\ \varepsilon\geq
0,$ (10c) $\displaystyle B^{\top}\lambda+D^{\top}\mu=f,\
e^{\top}\lambda+\eta=0,$ (10d) $\displaystyle\lambda,\mu,\eta\geq 0,\
\eqref{eq:linearized-compl}.$ (10e)

Exploiting Theorem 4 of [7], the existence of a finite and exact value for the
penalty parameter $\kappa$ is guaranteed. Moreover, Proposition 1 of [7]
ensures that the sets of optimal solutions of (10) and (8) are the same. In
addition, the polynomial size of $\kappa$ is shown by Theorem 22 of [9].
Finally, Problem (10) is the KKT reformulation of the bilevel problem (7). ∎

######  Corollary 2.3.

There is a polynomial-sized penalty parameter $\kappa>0$ so that the following
holds. For every bilevel optimal solution $(x,y)$ to Problem (1), the point
$(x,y,0)$ is bilevel optimal to Problem (7) with the same objective value. For
every bilevel optimal point $(x,y,\varepsilon)$ of Problem (7), the point
$(x,y)$ is bilevel optimal to Problem (1) with the same objective value.

######  Remark 2.4.

  1. (a)

Let us note that the reformulations of mixed-integer linear programs as linear
bilevel problems presented in [15] and [1] use two strategies for enforcing
that continuous lower-level variables are binary. One is based on coupling
constraints of the form $v=0$, where $v$ is an auxiliary lower-level variable,
and the other one corresponds to an exact penalization of these coupling
constraints and therefore requires no coupling constraints in the final model.
Our approach to penalize general coupling constraints in linear bilevel
problems uses the same idea in the final step.

  2. (b)

The derivations in this section suggest that one could potentially extend the
class of bilevel problems to which the results can be applied. First, the
results in [9] are valid for convex mixed-integer quadratic programs (MIQPs).
Hence, this would allow for convex-quadratic objective functions of the leader
and further integrality constraints in $X$. Moreover, we can also allow for
convex-quadratic (but still continuous) problems in the lower level since
their KKT conditions lead to polyhedral constraints in the KKT reformulation
(8). However, the big-$M$s from [4] cannot be used directly anymore. Given
that valid big-$M$s can also be found for quadratic programs in the lower
level, the most general class of bilevel problems to which our results could
be applied are those with convex MIQPs for the leader and convex QPs for the
follower.

  3. (c)

Due to the finiteness of the penalty parameter in Corollary 2.3, the original
bilevel problem with coupling constraints can also be solved by a finite
sequence of bilevel problems without coupling constraints if we follow
standard ideas of penalty methods.

##  3\. Discussion

It has been known for at least 25 years that coupling constraints in linear
bilevel problems can lead to disconnected feasible sets and that this is not
possible without coupling constraints. However, we prove that there is no
difference between these two types of linear bilevel problems on the level of
optimal solutions. While, on the one hand, this closes a gap in the literature
on linear bilevel optimization, it, on the other hand, also has some practical
consequences. Many theoretical statements in the literature on bilevel
optimization are made and shown for problems without coupling
constraints—either simply for the ease of presentation or due to a lack of a
proof for the case with coupling constraints. This note now allows for
carrying over some of these results by transforming the given problem having
coupling constraints into one without. We also point to the open question of
efficient computation of the penalty parameter required in the last section.
While we prove that it is polynomial-sized (in the bit-encoding length of the
data of the given problem), the question on how to compute it in polynomial
time is still open. Finally, the analogue question regarding the impact of
coupling constraints is, to the best of our knowledge, still open for
pessimistic bilevel problems.

## Acknowledgements

The second and third author acknowledge the support by the German
Bundesministerium für Bildung und Forschung within the project “RODES”
(Förderkennzeichen 05M22UTB).

## References

  * [1] Charles Audet, Pierre Hansen, Brigitte Jaumard and Gilles Savard  “Links between linear bilevel and mixed 0–1 programming problems”  In _Journal of Optimization Theory and Applications_ 93.2 Springer, 1997, pp. 273–300  DOI: 10.1023/A:1022645805569
  * [2] Omar Ben-Ayed and Charles E. Blair  “Computational Difficulties of Bilevel Linear Programming”  In _Operations Research_ 38.3 INFORMS, 1990, pp. 556–560  DOI: 10.1287/opre.38.3.556
  * [3] H. P. Benson  “On the structure and properties of a linear multilevel programming problem”  In _Journal of Optimization Theory and Applications_ 60.3, 1989, pp. 353–373  DOI: 10.1007/BF00940342
  * [4] Christoph Buchheim  “Bilevel linear optimization belongs to NP and admits polynomial-size KKT-based reformulations”  In _Operations Research Letters_ 51.6, 2023, pp. 618–622  DOI: 10.1016/j.orl.2023.10.006
  * [5] Stephan Dempe and Joydeep Dutta  “Is bilevel programming a special case of a mathematical program with complementarity constraints?”  In _Mathematical Programming_ 131.1-2 Springer, 2012, pp. 37–48  DOI: 10.1007/s10107-010-0342-1
  * [6] “Bilevel Optimization”, Springer Optimization and Its Applications 161  Cham: Springer, 2020  DOI: 10.1007/978-3-030-52119-6
  * [7] Mohammad Javad Feizollahi, Shabbir Ahmed and Andy Sun  “Exact augmented Lagrangian duality for mixed integer linear programming”  In _Mathematical Programming_ 161.1-2 Springer ScienceBusiness Media LLC, 2016, pp. 365–387  DOI: 10.1007/s10107-016-1012-8
  * [8] Michael R. Garey and David S. Johnson  “Computers and Intractability”, 1979
  * [9] Xiaoyi Gu, Shabbir Ahmed and Santanu S. Dey  “Exact Augmented Lagrangian Duality for Mixed Integer Quadratic Programming”  In _SIAM Journal on Optimization_ 30.1 Society for Industrial & Applied Mathematics (SIAM), 2020, pp. 781–797  DOI: 10.1137/19m1271695
  * [10] Pierre Hansen, Brigitte Jaumard and Gilles Savard  “New branch-and-bound rules for linear bilevel programming”  In _SIAM Journal on scientific and Statistical Computing_ 13.5 SIAM, 1992, pp. 1194–1217  DOI: 10.1137/0913069
  * [11] Robert G. Jeroslow  “The polynomial hierarchy and a simple model for competitive analysis”  In _Mathematical Programming_ 32.2 Springer, 1985, pp. 146–164  DOI: 10.1007/BF01586088
  * [12] Thomas Kleinert  “Algorithms for Mixed-Integer Bilevel Problems with Convex Followers”, 2021  URL: https://opus4.kobv.de/opus4-trr154/frontdoor/index/index/docId/383
  * [13] Patrice Marcotte and Gilles Savard  “Bilevel Programming: A Combinatorial Perspective”  In _Graph Theory and Combinatorial Optimization_ Boston, MA: Springer US, 2005, pp. 191–217  DOI: 10.1007/0-387-25592-3_7
  * [14] Luis Vicente, Gilles Savard and Joaquim Júdice  “Descent approaches for quadratic bilevel programming”  In _Journal of Optimization Theory and Applications_ 81.2 Springer, 1994, pp. 379–399  DOI: 10.1007/BF02191670
  * [15] Luis Vicente, Gilles Savard and Joaquim Júdice  “Discrete linear bilevel programming problem”  In _Journal of Optimization Theory and Applications_ 89.3, 1996, pp. 597–614  DOI: 10.1007/BF02275351


# Dynamical excitations of one-dimensional Fulde-Ferrell pairing Fermi
superfluid

Peng Zou1 Huaisong Zhao1<EMAIL_ADDRESS>Feng Yuan1 Shi-Guo Peng2
<EMAIL_ADDRESS>1College of Physics, Qingdao University, Qingdao 266071,
China  2State Key Laboratory of Magnetic Resonance and Atomic and Molecular
Physics, Innovation Academy for Precision Measurement Science and Technology,
Chinese Academy of Sciences, Wuhan 430071, China

###### Abstract

We theoretically investigate a one-dimensional Fulde-Ferrell Fermi superfluid
at a finite effective Zeeman field $h$, and study entire dynamical excitations
related to density perturbation. By calculating the density dynamic structure
factor, we find anisotropic dynamical excitations in both collective modes and
single-particle excitations. Along the direction of centre-of-mass momentum
$p$, there are two obvious gapless collective modes with different speed. The
lower collective modes is from the usual gauge symmetry breaking and has a
larger speed than the one in the negative direction of $p$. The higher one is
due to the direction spontaneous symmetry breaking of centre-of-mass momentum
$p$, and separates two kinds of single-particle excitations in the positive
$p$ direction. However, this higher mode disappears in the opposite direction
of $p$, where two single-particle excitations overlap with each other. These
signals of dynamical excitations can do help to distinguish Fulde-Ferrell
superfluid from the conventional Bardeen-Cooper-Schrieffer superfluid in the
future experiment.

##  I Introduction

Since the discovery of superconducting phenomena in mercury by K. Onnes in
1911, it is realized that particles with opposite momentum and spin can
generate a molecular Cooper pair carrying zero center-of-mass (COM) momentum
to decrease energy. This conventional Bardeen-Cooper-Schrieffer (BCS)
superfluidity originates from the Bose-Einstein-Condensation (BEC) of Cooper
pairs at zero momentum. Later it is theoretically predicted by Fulde and
Ferrel (FuldePR1964, ) and Larkin and Ovchinnikov (Larkin1964, ) that Cooper
pairs can also carry a finite COM momentum, and condense at a nonzero momentum
$p$ with order parameter in either a plane-wave FF-type
$\Delta\left(\bm{r}\right)=\Delta e^{i\bm{p}\centerdot\bm{r}}$ or a standing-
wave LO-type $\Delta\left(\bm{r}\right)=\Delta{\rm
cos}\left(\bm{p}\centerdot\bm{r}\right)$. In such an ansatz, the two
mismatched Fermi surfaces of different spin components can overlap, thereby
supporting a spatially inhomogeneous superfluidity. These pairing states
generate an exotic Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) superfluidity
(CasalbuoniRMP2004, ).

These FFLO superfluids have been actively searched for almost six decades. In
condensed matter physics, some strong signals of FFLO states come from the
research in heavy fermions (RadovanNAT2003, ; GloosPRL1993, ; HuxleyCM1993, ;
BianchiPRL2003, ; MartinPRB2005, ) and organic superconductors
(Croitoru2017CM, ), which are only indirect experimental evidences. Recently
several experimental works reported evidence of pair-density-wave in high-Tc
superconductor by scanning tunnelling microscopy technique (LiuNAT2023, ;
ZhaoNAT2023, ; AishwaryaNAT2023, ), which ignites the passion again and
provides new possibilities to search and study FFLO state. An ultracold atomic
Fermi gas has proven to be an ideal tabletop system for the pursuit of FFLO
superfluidity (RadzihovskyRPP2010, ). The FFLO state is thought to be very
fragile in three dimensions and has quite narrow parameter space. Both lower
spatial dimension (LiuPRA2007, ) and spin-orbit coupling effect (HuNJP2013, ;
LiuPRA2013, ) in Fermi superfluid have been reported that they can expand the
parameter space and theoretically investigate the possible phase diagrams
related to FFLO state. Also a theoretical strategy via a Dark-State Control of
Feshbach resonance was proposed to realize FF superfluid in ultracold atomic
gases (HePRL2018, ).

In fact it is expected that the full dynamical excitations of a certain matter
state can be utilized to check the existence of this state, or distinguish it
from other states. Dynamic structure factor is an important many-body physical
quantity, and includes rich information related to dynamical excitations of
the system (Bloch2008rmp, ; Giorgini2008rmp, ). Experimentally dynamic
structure factor can be directly measured by a two-photon Bragg scattering
technique, which had been used to investigate dynamical excitations of the
BCS-BEC crossover Fermi superfluid, including the single-particle excitations
(VeeravalliPRL2008, ), collective Goldstone phonon mode (HoinkNP2017, ),
second sound (LiSCI2022, ) and Higgs mode (DykeArxiv2023, ). A spin-charge
separation of repulsive one-dimensional (1D) Fermi gas has also been studied
by this technique (SenaratneSCI2022, ). Thus it is interesting to study the
dynamical excitations of an FFLO Fermi superfluid, and find its dynamical
characters related to its symmetry structure.

In this paper, we theoretically investigate a 1D spin-polarized Fermi
superfluid in an FF-type pairing state, and discuss its entire dynamical
excitations by numerically calculating the dynamic structure factor of this
system with random phase approximation (AndersonPR1958, ; LiuPRA2004, ). We
will initially discuss the state of equation of the system in FF superfluid,
and then introduce all possible dynamical excitations, which may provide some
background knowledge to anisotropic Josephson effect related to FFLO state.
This paper is organized as follows. In the next section, we use the motion
equation of Green’s function to introduce the microscopic model of a 1D spin-
polarized interacting Fermi gas, and outline the mean-field approximation, and
then show how to calculate the response function with random phase
approximation. We give results of the dynamic structure factor of FF
superfluids in Sec. III, and we give our conclusions and outlook in Sec. IV.
Some calculation details are listed in the Appendix.

##  II Methods

###  II.1 Model and Hamiltonian

We consider a uniform 1D spin-polarized Fermi gases with s-wave contact
interaction. The system can be described by a model Hamiltonian

$H=\sum_{k\sigma}\left(\epsilon_{k}-\mu_{\sigma}\right)c_{k\sigma}^{\dagger}c_{k\sigma}+g_{1D}\sum_{pkk^{\prime}}c_{k\uparrow}^{\dagger}c_{p-k\downarrow}^{\dagger}c_{p-k^{\prime}\downarrow}c_{k^{\prime}\uparrow},$
(1)

here $\epsilon_{k}=k^{2}/2m$ is the dispersion relation of spin-$\sigma$ free
particles with mass $m$ in reference to the chemical potential $\mu_{\sigma}$,
and $c_{k\sigma}(c_{k\sigma}^{\dagger})$ is their annihilation (generation)
operator in momentum representation. Here and in the following we have set
physical quantities $\hbar=k_{B}=1$ for convenience. The interaction
$g_{1D}=-\gamma n_{0}/m$ describes an attractive interplay between opposite
spin components with a dimensionless quantity $\gamma$. Since we consider a
uniform system with bulk density $n_{0}$, the inverse of Fermi wave vector
$k_{F}=\pi n_{0}$/2 and Fermi energy $E_{F}=k_{F}^{2}/2m$ are used as length
and energy units, respectively. Usually the difference of chemical potential
is used to define an effective Zeeman field
$h=\left(\mu_{\uparrow}-\mu_{\downarrow}\right)/2$, and
$\mu=\left(\mu_{\uparrow}+\mu_{\downarrow}\right)/2$ is the average chemical
potential.

At zero temperature $T=0$, usually the system will comes into a superfluid
state, in which two opposite-spin atoms form a molecular Cooper pair to
decrease energy. Interestingly, an FF-type superfluid, who has a non-zero COM
momentum in Cooper pair, possibly turns out and even becomes the ground state
at a non-zero effective Zeeman field $h$. In a standard mean-field treatment,
the order parameter of FF superfluid can be expressed in a plane-wave form
$\Delta\left(x\right)\equiv\Delta e^{ipx}$, where COM momentum $p$ and pairing
gap $\Delta=g_{1D}\sum_{{}_{k}}\left\langle
c_{-k+p/2\downarrow}c_{k+p/2\uparrow}\right\rangle$ are two important degrees
of freedom of order parameter $\Delta\left(x\right)$. Within this
approximation, a mean-field Hamiltonian is expressed as

$\begin{array}[]{cl}H_{\rm{mf}}&=\sum_{k}\left(\xi_{k}-h\sigma_{z}\right)c_{k\sigma}^{\dagger}c_{k\sigma}-\Delta^{2}/g_{1D}\\\
&-\sum_{k}\Delta\left(c_{-k+p/2\downarrow}c_{k+p/2\uparrow}+h.c.\right),\end{array}$
(2)

where $\xi_{k}=\epsilon_{k}-\mu$. The exact solution of mean-field Hamiltonian
$H_{{\rm mf}}$ can be obtained by motion equation of Green’s functions.
Finally we get three Green’s functions, whose expressions are listed below

$G_{1}\equiv\left\langle\left\langle
c_{k+p/2\uparrow}|c_{k+p/2\uparrow}^{\dagger}\right\rangle\right\rangle=\underset{l}{\sum}\frac{\left[G_{1}\right]_{k}^{l}}{\omega-
E_{k}^{l}},$ (3)

$G_{2R}\equiv\left\langle\left\langle
c_{-k+p/2\downarrow}^{\dagger}|c_{-k+p/2\downarrow}\right\rangle\right\rangle=\underset{l}{\sum}\frac{\left[G_{2R}\right]_{k}^{l}}{\omega-
E_{k}^{l}},$ (4)

$\varGamma\equiv\left\langle\left\langle
c_{k+p/2\uparrow}|c_{-k+p/2\downarrow}\right\rangle\right\rangle=\underset{l}{\sum}\frac{\left[\varGamma\right]_{k}^{l}}{\omega-
E_{k}^{l}},$ (5)

The double brackets in above equations are used to define Fourier
transformation of double-time Green’s function. Here it should be emphasized
that the expression of spin-down Green’s functions $G_{2R}$ is different from
the one of spin-up Green’s function $G_{1}$. Their definitions closely depend
on their way coupled to pairing Green’s function $\Gamma$ when solving motion
equation of Green’s function. All expressions related
to$\left[G_{1}\right]_{k}^{l}$, $\left[G_{2R}\right]_{k}^{l}$ and
$\left[\Gamma\right]_{k}^{l}$ will be listed in the last appendix. $l=1,2$
denotes two branches of quasi-particle energy spectrum $E_{k}^{(1)}$ and
$E_{k}^{(2)}$,

$E_{k}^{(1,2)}=\frac{kp}{2m}-h\pm E_{k},$ (6)

where $E_{k}=\sqrt{\xi_{kp}^{2}+\Delta^{2}}$ and
$\xi_{kp}=\epsilon_{k}+\epsilon_{p}/4-\mu$. In typical parameters used in this
paper, we find that the value of $E_{k}^{(2)}$ is always negative, while
$E_{k}^{\left(1\right)}$ can be either positive or negative. The distributions
of these two spectra are shown in Fig. 1, in which $E_{k}^{\left(1\right)}$ is
negative when $-2.17k_{F}<k<0.12k_{F}$.

Figure 1: Two branches of single-particle excitation spectrum at interaction
strength $\gamma=3$ and effective Zeeman field $h=1.2E_{F}$. Here COM momentum
$p=1.18k_{F}$.

The mean-field thermodynamic potential of the system reads

$\begin{array}[]{cl}\Omega=&-\Delta^{2}/g_{1D}+\sum_{k}\left(\xi_{k}-E_{k}\right)\\\
&+\sum_{k}T\left[\rm{ln}\it{f}\left(-E_{k}^{(\rm{1})}\right)+\rm{ln}\it{f}\left(E_{k}^{(\rm{2})}\right)\right],\end{array}$
(7)

which is connected to the free energy by $F=\Omega+\mu N$. Here
$f\left(x\right)=1/\left(e^{x/T}+1\right)$ is the Fermi-Dirac distribution
function at temperature $T$. In this paper, we focus our discussion on an
almost zero temperature ($T=0.01T_{F}$), to avoid an unnecessary numerical
divergence induced by zeros of $E_{k}^{(1)}$.

###  II.2 State of equations

The values for average chemical potential $\mu$, amplitude of order parameter
$\Delta$ and COM momentum $p$ can be respectively determined with minimization
of thermodynamic potential $\Omega$ in Eq. 7 to them, namely
$N=-\partial\Omega/\partial\mu$, $\partial\Omega/\partial\Delta=0$ and
$\partial\Omega/\partial p=0$. These three relations respectively give the
total particle number equation

$N=\sum_{k}\left(1-\frac{\xi_{kp}}{E_{k}}\right)+\sum_{k}\frac{\xi_{kp}}{E_{k}}\left[f\left(E_{k}^{(1)}\right)+f\left(-E_{k}^{(2)}\right)\right],$
(8)

pairing gap equation

$\frac{\Delta}{g_{1D}}=\sum_{k}\frac{\Delta}{2E_{k}}\left[f\left(E_{k}^{(1)}\right)+f\left(-E_{k}^{(2)}\right)-1\right],$
(9)

and COM momentum equation

$\sum_{k}\left[p\left(1-\frac{\xi_{kp}}{E_{k}}\right)+\left(2k+\frac{\xi_{kp}}{E_{k}}p\right)f\left(E_{k}^{(1)}\right)-\left(2k-\frac{\xi_{kp}}{E_{k}}p\right)f\left(-E_{k}^{(2)}\right)\right]=0.$
(10)

The value of $\mu$, $\Delta$ and $p$ should be self-consistently solved with
Eqs. 8, 9 and 10.

Figure 2: State of equation at interaction strength $\gamma=3$. Three panels
show the curves of (a) free energy, (b) COM momentum and (c) average chemical
potential and order parameter at different effective Zeeman field,
respectively.  Figure 3: The critical effective Zeeman field $h_{c}$ at
different interaction strength $\gamma$.

The possible stable state of a system is determined by all minima in free
energy $F$. Generally there are three possible states here. Besides the
trivial normal state, who always has a zero pairing gap $\Delta$ and the
largest free energy, the conventional BCS superfluid and FF superfluid are the
other two states with relatively lower free energy. As shown in Fig. 2, at a
small $h$, the conventional BCS superfluid is the ground state of the system,
whose COM momentum $p=0$. The value of $h$ can hardly influence the chemical
potential $\mu$ and pairing gap $\Delta$, which can be thought as the analogy
of the Meissner effect in the superfluid. Interestingly an FF superfluid with
a non-zero COM momentum $p$ begins to turn out when $h$ is large enough, whose
free energy is larger than the one of BCS superfluid. When $h$ comes to a
critical value $h_{c}$ (a little smaller than $1.0E_{F}$ here), this FF
superfluid shares the same free energy as BCS superfluid. Further increasing
$h$, FF superfluid will replace BCS superfluid to be the new ground state of
the system for $h>h_{c}$. This phase transition had been introduced in Ref.
(LiuPRA2007, ). Here $\Delta$ and $p$ are two necessary degree of freedom of
in a FF-type order parameter $\Delta\left(x\right)=\Delta e^{ipx}$. Pairing
gap $\Delta$ is related to the conventional gauge symmetry breaking between
normal state and superfluidity, and the direction of COM momentum $p$ is
related to another spontaneous symmetry breaking since $\pm p$ corresponding
the same free energy. While $\Delta$ always decreasing with $h$, $p$ shows a
monotonically increasing behaviour with $h$. The different dependence behavior
of $\Delta$ and $p$ with effective Zeeman field $h$ goes on reflecting that
they come from different symmetry breaking and may induce different collective
modes. It is easy to know that the phase transition shown in Fig. 2 is a
first-order one, which is manifested by the discontinuous behavior of $\mu$,
$\Delta$ and $p$. We have checked that the same phase transition also happens
at different interaction strength $\gamma$ with different critical Zeeman
field $h_{c}$, whose value is shown in Fig. 3. A larger interaction strength
$\gamma$ requires a bigger critical effective Zeeman field $h_{c}$ to make the
system come into the FF superfluid state.

###  II.3 Calculation of dynamical excitations

When an interacting system comes into a superfluid state, usually there will
be four typical densities. Besides the normal spin-up density
$n_{1}=\left\langle\psi_{\uparrow}^{\dagger}\psi_{\uparrow}\right\rangle$ and
spin-down density
$n_{2}=\left\langle\psi_{\downarrow}^{\dagger}\psi_{\downarrow}\right\rangle$,
the pairing physics of opposite-spin atoms generates the other anomalous
density $n_{3}=\left\langle\psi_{\downarrow}\psi_{\uparrow}\right\rangle$ and
its conjugate counterpart
$n_{4}=\left\langle\psi_{\uparrow}^{\dagger}\psi_{\downarrow}^{\dagger}\right\rangle$.
These pairs form phase coherent Cooper pairs with zero or finite COM momentum
$p$. The interaction between particles makes these four densities couple
closely with each other. Any fluctuation in each kind of density will
influence other densities and generate an obvious density fluctuation of them.
Also any weak perturbation potential $V_{\textrm{pert}}$ will generate density
fluctuations $\delta n$, and they are connected with each other by response
function $\chi$, namely $\delta n=\chi V_{\textrm{pert}}$, in the frame of
linear response theory.

The entire dynamical excitations of the system consist of possible collective
excitations and single-particle excitations, which are also mainly connected
to the physical properties of Cooper pairs. The phase spontaneously breaking
of them generates a gapless Goldstone collective mode, while the breaking of
Cooper pairs forming parts of single-particle excitation. These dynamical
excitations can be well described by the density dynamic structure factor,
which is from the imaginary part of related response functions $\chi$. The
direct calculation of dynamic structure factor suffers from the many-body
difficulty. A potential approximation to overcome this problem is the random
phase approximation, which is firstly brought by P. W. Anderson and has be
verified to be a qualitatively reliable way to study dynamical excitations. In
our previous work, we had introduced how to calculate dynamic structure factor
with random phase approximation (GaoPRA2023, ; ZhaoPRA2023, ). This
approximation finds the connection between the beyond mean-field response
function $\chi$ and its mean-field approximation $\chi^{0}$, whose calculation
is relatively easier, by the following equation

$\chi=\frac{\chi^{0}}{1-\chi^{0}M_{I}g_{1D}},$ (11)

where the constant matrix $M_{I}=\sigma_{0}\otimes\sigma_{x}$ is the direct
product of unit matrix $\sigma_{0}$ and Pauli matrix $\sigma_{x}$, reflecting
the coupling situation of four types of densities. $\chi^{0}$ is a four
dimension matrix in mean-field level, and its $ij$ matrix element
$\chi_{ij}^{0}$ reflects the interaction-induced coupling between density
$n_{i}$ and $n_{j}$. When the system comes into the FF superfluid, $\chi^{0}$
has nine independent matrix elements, namely

$\chi^{0}=\left[\begin{array}[]{cccc}\chi_{11}^{0}&\chi_{12}^{0}&\chi_{13}^{0}&\chi_{14}^{0}\\\
\chi_{12}^{0}&\chi_{22}^{0}&\chi_{23}^{0}&\chi_{24}^{0}\\\
\chi_{14}^{0}&\chi_{24}^{0}&-\chi_{12}^{0}&\chi_{34}^{0}\\\
\chi_{13}^{0}&\chi_{23}^{0}&\chi_{43}^{0}&-\chi_{12}^{0}\end{array}\right].$
(12)

After Fourier transformation to above response function, we obtain the
expression of all matrix elements in the momentum-energy representation.

With Eqs. 11 and 12, we can obtain the expression of total-density response
function $\chi_{n}\equiv\chi_{11}+\chi_{22}+\chi_{12}+\chi_{21}$. Based on the
fluctuation and dissipation theorem, its imaginal part is connected to density
dynamical structure factor by

$S_{n}=-\frac{1}{\pi}\frac{1}{1-e^{-\omega/T}}{\rm Im}\left(\chi_{n}\right).$
(13)

##  III Results

Figure 4: The density dynamic structure factor $S_{n}(q,\omega)$ of FF
superfluid at interacting strength $\gamma=3$. (L) the lower collective phonon
mode along the COM momentum $p$ direction, (H) the higher collective phonon
mode, (N) the collective phonon mode in the opposite direction of $p$.

We numerically calculate the density dynamic structure factor
$S_{n}(q,\omega)$ of an FF superfluid at interaction strength $\gamma=3$. We
define the direction of COM momentum $p$ is along the positive direction of
transferred momentum $q$. The results are shown in Fig. 4. Intuitively we see
an anisotropic dynamical behavior. The dynamical excitation at a positive
transferred momentum $q$ is different from the one at a negative $q$. This is
due to the direction dependence of the non-zero COM momentum $p$ in FF
superfluid. Along the direction of COM momentum $p$ (namely $q>0$), we see two
kinds of gapless collective mode (curves marked by (L) and (H)), and two-
separated regimes of single-particle excitation. When $q<0$, the higher
gapless collective mode disappears, and two kinds of single-particle
excitations overlap with each other. The speed of gapless collective phonon
mode (curve marked by (N)) is a little smaller than that in the positive
direction. In the following, we will separately introduce collection modes and
single-particle excitations.

###  III.1 Collective excitation

The origin of collective mode is closely related to the symmetry breaking of a
certain matter state. A gapless collective excitation comes from a certain
spontaneous symmetry breaking of the system. It is interesting to notice that
there are two gapless collective modes at a positive transferred momentum $q$,
which are also displayed at Fig. 4. The lower mode is the conventional
collective Goldstone phonon mode (marked by (L)), which requires the lowest
excitation energy among all possible excitations. Its physical origin is due
to the gauge symmetry breaking of the phase of order parameter $\Delta(x)$.
The higher gapless collective mode (marked by (H)) is due to the symmetry
breaking of the direction of COM momentum $p$, which is continuous symmetry
breaking in higher dimension system but not continuous in 1D system due to its
specific spatial dimension. A similar mode as the higher collective mode here
is also reported in a LO-type superconductor with order parameter
$\Delta\left(x\right)=\Delta\rm{cos}\left(\it{px}\right)$, which is called the
gapless Higgs mode since its amplitude of order parameter displays spatial
periodic variation. However, this should be different from the one in FF
superfluid, whose amplitude is always a constant value to keep continuous
translational symmetry. So the higher collective mode in FF superfluid we
argue that it is better to be called a gapless phonon-like mode, instead of a
gapless Higgs mode.

Figure 5: Speeds of three collective modes at different effective Zeeman field
$h$. $v_{L}$ and $v_{H}$ are respectively the speed of lower and higher branch
collective mode along the direction of COM momentum $p$, while $v_{N}$ is the
speed of collective mode in the opposite direction of $p$.

The other difference from the BCS superfluid is that the lowest collective
phonon mode displays an anisotropic excitation between the positive and
negative direction, which is due to the direction dependence of COM momentum
$p$ in FF superfluid. The similar anisotropic dynamics of phonon mode is also
reported in a spin-polarized Fermi system in square optical lattice
(HeikkinenPRA2011, ). From Fig. 5, it is easy to know that the velocity of
phonon mode in a negative $p$ direction (marked by (N)) is always smaller than
the one in the positive $q$ direction. Also the absolute value of phonon
velocity decreases with $h$ in the negative $q$ direction.

The effective Zeeman field $h$-dependence of the speed for three collective
modes is shown in Fig. 5, which displays that both $v_{L}$ and $v_{N}$
decrease with $h$ and slowly goes to zero at a large enough $h$. Finally
$\Delta$ also touch zero and the superfluid disappears when $h$ is around
$2.0E_{F}$. However, $v_{H}$ always rises with $h$, and this behavior is
consistent with that of COM momentum $p$ (panel (b) in Fig. 2). The similar
$h$-dependence between $v_{H}$ and $p$ again demonstrates a close connection
between them. Also the different $h$-dependence between $v_{L}$ and $v_{N}$
demonstrates that the two gapless collective modes come from different
symmetry breaking mechanism, and are related to different fluctuation of order
parameter.

###  III.2 Single-particle excitation

Now we discuss the pair-breaking excitation, which is an important part of
single-particle excitation and takes up large regimes in Fig. 4. This
excitation is usually a continuous excitation and its lowest excitation energy
is determined by the single-particle spectrum (Eq. 6). To understand all
possible ways of pair-breaking excitation, it is better to understand this
physics by part of expression in response function $\chi^{0}$, namely

$\frac{f\left(E_{p}^{l}\right)-f\left(E_{p+q}^{l^{\prime}}\right)}{i\omega_{n}+E_{p}^{l}-E_{p+q}^{l^{\prime}}}.$
(14)

Here $E_{p}^{l}$ should consider all possible combinations of single-particle
spectrum. As shown in Fig. 1, $E_{k}^{(2)}$ is always negative, and
$E_{k}^{\left(1\right)}$ is positive except in a narrow regime
$-2.17k_{F}<k<0.12k_{F}$. The possible pair-breaking excitation may happen
requiring that the numerator of Eq. 14 can not be zero, and also
$E_{p}^{l}-E_{p+q}^{l^{\prime}}<0$ to ensure a positive excitation energy.

Finally we find two possible pair-breaking excitations, namely
$\left|E_{k}^{(1)}-E_{k+q}^{(1)}\right|$ (11-type excitation) at a limited
regime of momentum $k$ and $\left|E_{k}^{(2)}-E_{k+q}^{(1)}\right|$(21-type
excitation) at a full regime of $k$. And their minima is the lowest energy to
break an FF-type Cooper pair, which are labeled by both pink and yellow dotted
lines in Fig. 4, respectively. It should emphasized that the 11-type
excitations (dotted yellow line) is a gapless pair-breaking excitation and is
absent in the conventional BCS superfluid. Honestly there is slightly bias
between predictions from yellow dotted lines and random phase approximation’s
prediction when $q$ is around $\pm 2.5k_{F}$, which may due to the limited
excitation regime of momentum $k$ in 11-excitation. The minimum of 21-type
excitation is labeled with a dotted pink line. For a positive transferred
momentum $q$, these two pair-breaking excitations are separated with each
other by just the higher gapless collective mode. However, they are mixed with
each other in a negative transferred momentum $q$. These differences can do
help to find the higher collective mode in future experiment, and distinguish
FF-type superfluid from BCS superfluid.

##  IV Conclusions and outlook

In summary, we theoretically calculate the dynamic structure factor of a 1D FF
superfluid with random phase approximation to study dynamical excitations of
the system. We find an anisotropic dynamical excitation between positive and
negative directions in both collective modes and single-particle excitations.
In the positive direction, we find two gapless collective modes. The lower one
comes from the spontaneous breaking of gauge symmetry, while the higher one
comes from the direction symmetry breaking of COM momentum $p$. The sound
speed in the positive direction is larger than the one in the negative
direction. There are two types of pair-breaking excitations, and one of them
is a gapless excitation, which is absent in the BCS superfluid. In the
positive direction, these two kinds of pair-breaking excitations are just
separated by the higher gapless collective mode, but overlap with each other
in the negative direction. These dynamical excitations can do help to
distinguish FF superfluid from the conventional BCS superfluid in future
experiment.

##  V Acknowledgements

We are grateful for discussions with Wei Yi, Xiaoquan Yu and Hui Hu. This
research was supported by National Natural Science Foundation of China under
Grants No. 11804177 (P.Z.), No. 11547034 (H.Z.), Grant No. 11974384 (S.-G.P.)
and 12374250 (S.-G.P.), and the National Key R&D Program under Grant No.
2022YFA1404102 (S.-G.P.).

##  VI Appendix

In this appendix, we list expressions of three Green’s functions and mean-
field response function $\chi^{0}$. The first Green’s functions is
$G_{1}\left(k,\omega\right)=\sum_{l}\left[G_{1}\right]_{k}^{l}/\left(\omega-
E_{k}^{l}\right)$, with

$\begin{array}[]{cc}\left[G_{1}\right]_{k}^{(1)}=\frac{E_{k}^{(1)}+\xi_{k-p/2}+h}{E_{k}^{(1)}-E_{k}^{(2)}},&\left[G_{1}\right]_{k}^{(2)}=-\frac{E_{k}^{(2)}+\xi_{k-p/2}+h}{E_{k}^{(1)}-E_{k}^{(2)}}.\end{array}$

The second one is
$G_{2R}\left(k,\omega\right)=\sum_{l}\left[G_{2R}\right]_{k}^{l}/\left(\omega-
E_{k}^{l}\right),$ with

$\begin{array}[]{cc}\left[G_{2R}\right]_{k}^{(1)}=\frac{E_{k}^{(1)}-\xi_{k+p/2}+h}{E_{k}^{(1)}-E_{k}^{(2)}},&\left[G_{2R}\right]_{k}^{(2)}=-\frac{E_{k}^{(2)}-\xi_{k+p/2}+h}{E_{k}^{(1)}-E_{k}^{(2)}}.\end{array}$

The third Green’s function is
$\varGamma\left(k,\omega\right)=\sum_{l}\left[\varGamma\right]_{k}^{l}/\left(\omega-
E_{k}^{l}\right),$ with

$\begin{array}[]{cc}\left[\Gamma\right]_{k}^{(1)}=-\frac{\Delta}{E_{k}^{(1)}-E_{k}^{(2)}},&\left[\Gamma\right]_{k}^{(2)}=\frac{\Delta}{E_{k}^{(1)}-E_{k}^{(2)}}.\end{array}$

The expressions of all nine independent matrix elements in mean-field response
function $\chi^{0}$ are respectively

$\chi_{11}^{0}=+\underset{pll^{\prime}}{\sum}\left[G_{1}\right]_{p}^{l}\left[G_{1}\right]_{p+q}^{l^{\prime}}\frac{f\left(E_{p}^{l}\right)-f\left(E_{p+q}^{l^{\prime}}\right)}{i\omega_{n}+E_{p}^{l}-E_{p+q}^{l^{\prime}}},$

$\chi_{12}^{0}=-\underset{pll^{\prime}}{\sum}\left[\varGamma\right]_{p}^{l}\left[\varGamma\right]_{p+q}^{l^{\prime}}\frac{f\left(E_{p}^{l}\right)-f\left(E_{p+q}^{l^{\prime}}\right)}{i\omega_{n}+E_{p}^{l}-E_{p+q}^{l^{\prime}}},$

$\chi_{13}^{0}=+\underset{pll^{\prime}}{\sum}\left[G_{1}\right]_{p}^{l}\left[\varGamma\right]_{p+q}^{l^{\prime}}\frac{f\left(E_{p}^{l}\right)-f\left(E_{p+q}^{l^{\prime}}\right)}{i\omega_{n}+E_{p}^{l}-E_{p+q}^{l^{\prime}}},$

$\chi_{14}^{0}=+\underset{pll^{\prime}}{\sum}\left[\varGamma\right]_{p}^{l}\left[G_{1}\right]_{p+q}^{l^{\prime}}\frac{f\left(E_{p}^{l}\right)-f\left(E_{p+q}^{l^{\prime}}\right)}{i\omega_{n}+E_{p}^{l}-E_{p+q}^{l^{\prime}}},$

$\chi_{22}^{0}=+\underset{kll^{\prime}}{\sum}\left[G_{2R}\right]_{k}^{l}\left[G_{2R}\right]_{k+q}^{l^{\prime}}\frac{f\left(E_{k}^{l}\right)-f\left(E_{k+q}^{l^{\prime}}\right)}{i\omega_{n}+E_{k}^{l}-E_{k+q}^{l^{\prime}}},$

$\chi_{23}^{0}=-\underset{kll^{\prime}}{\sum}\left[\varGamma\right]_{k}^{l}\left[G_{2R}\right]_{k+q}^{l^{\prime}}\frac{f\left(E_{k}^{l}\right)-f\left(E_{k+q}^{l^{\prime}}\right)}{i\omega_{n}+E_{k}^{l}-E_{k+q}^{l^{\prime}}},$

$\chi_{24}^{0}=-\underset{kll^{\prime}}{\sum}\left[G_{2R}\right]_{k}^{l}\left[\varGamma\right]_{k+q}^{l^{\prime}}\frac{f\left(E_{k}^{l}\right)-f\left(E_{k+q}^{l^{\prime}}\right)}{i\omega_{n}+E_{k}^{l}-E_{k+q}^{l^{\prime}}},$

$\chi_{34}^{0}=+\underset{kll^{\prime}}{\sum}\left[G_{2R}\right]_{k}^{l}\left[G_{1}\right]_{k+q}^{l^{\prime}}\frac{f\left(E_{k}^{l}\right)-f\left(E_{k+q}^{l^{\prime}}\right)}{i\omega_{n}+E_{k}^{l}-E_{k+q}^{l^{\prime}}},$

$\chi_{43}^{0}=+\underset{kll^{\prime}}{\sum}\left[G_{1}\right]_{k}^{l}\left[G_{2R}\right]_{k+q}^{l^{\prime}}\frac{f\left(E_{k}^{l}\right)-f\left(E_{k+q}^{l^{\prime}}\right)}{i\omega_{n}+E_{k}^{l}-E_{k+q}^{l^{\prime}}}.$

## References

  * (1) P. Fulde and R. A. Ferrell, Superconductivity in a Strong Spin-Exchange Field, Phys. Rev. 135, A550 (1964).
  * (2) A. I. Larkin and Y. N. Ovchinnikov, Inhomogeneous state of superconductors, Zh. Eksp. Teor. Fiz. 47, 1136 (1964); [Sov. Phys. JETP 20, 762 (1965)].
  * (3) For reviews on the FFLO states, see, for examples, R. Casalbuoni and G. Nardulli, Inhomogeneous superconductivity in condensed matter and QCD, Rev. Mod. Phys. 76, 263 (2004); J. J. Kinnunen, J. E. Baarsma, J.-P. Martikainen and P. Törmä, The Fulde–Ferrell–Larkin–Ovchinnikov state for ultracold fermions in lattice and harmonic potentials: a review, Rep. Prog. Phys. 81, 046401 (2018).
  * (4) Y. Matsuda and H. Shimahara, Fulde–Ferrell–Larkin–Ovchinnikov State in Heavy Fermion Superconductors, J. Phys. Soc. Jpn. 76, 051005 (2007).
  * (5) H. A. Radovan, N. A. Fortune, T. P. Murphy, S. T. Hannahs, E. C. Palm, S. W. Tozer, and D. Hall, Magnetic enhancement of superconductivity from electron spin domains, Nature 425, 51 (2003).
  * (6) K. Gloos, R. Modler, H. Schimanski, C. D. Bredl, C. Geibel, F. Steglich, A. I. Buzdin, N. Sato, and T. Komatsubara, Possible formation of a nonuniform superconducting state in the heavy-fermion compound $UPd_{2}Al_{3}$, Phys. Rev. Lett. 70, 501 (1993).
  * (7) A. D. Huxley, C. Paulson, O. Laborde, J. L. Tholence, D. Sanchez, A. Junod, and R. Calemczuk, Flux pinning, specific heat and magnetic properties of the laves phase superconductor $CeRu_{2}$, J. Phys.: Condens. Matter 5, 7709 (1993).
  * (8) A. Bianchi, R. Movshovich, C. Capan, P. G. Pagliuso, and J. L. Sarrao, Possible Fulde-Ferrell-Larkin-Ovchinnikov Superconducting State in $CeCoIn_{5}$, Phys. Rev. Lett. 91, 187004 (2003).
  * (9) C. Martin, C. C. Agosta, S. W. Tozer, H. A. Radovan, E. C. Palm, T. P. Murphy, and J. L. Sarrao, Evidence for the Fulde-Ferrell-Larkin-Ovchinnikov state in $CeCoIn_{5}$ from penetration depth measurements, Phys. Rev. B 71, 020503(R) (2005).
  * (10) M. Croitoru and A. Buzdin, In search of unambiguous evidence of the Fulde-Ferrell-Larkin-Ovchinnikov state in quasi-low-dimensional superconductors, Condens. Matter 2, 30 (2017).
  * (11) Y. Liu, T. Wei, G. He, Y. Zhang, Z. Wang & J. Wang, Pair density wave state in a monolayer high-Tc iron-based superconductor, Nature 618, 934 (2023).
  * (12) H. Zhao, R. Blackwell, M. Thinel, T. Handa, S. Ishida, X. Zhu, A. Iyo, H. Eisaki, A. N. Pasupathy & K. Fujita, Smectic pair-density-wave order in $EuRbFe_{4}As_{4}$, Nature 618, 940 (2023).
  * (13) A. Aishwarya, J. M.-Mann, A. Raghavan, L. Nie, M. Romanelli, S. Ran, S. R. Saha, J. Paglione, N. P. Butch, E. Fradkin & V. Madhavan, Magnetic-field-sensitive charge density waves in the superconductor $UTe_{2}$, Nature 618, 928 (2023).
  * (14) L. Radzihovsky and D. E. Sheehy, Imbalanced Feshbach-resonant Fermi gases, Rep. Prog. Phys. 73, 076501 (2010).
  * (15) X.-J. Liu, H. Hu, and P. D. Drummond, Fulde-Ferrell-Larkin-Ovchinnikov states in one-dimensional spin-polarized ultracold atomic Fermi gases, Phys. Rev. A 76, 043605 (2007).
  * (16) H. Hu and X.-J. Liu, Fulde–Ferrell superfluidity in ultracold Fermi gases with Rashba spin–orbit coupling, New J. Phys. 15, 093037 (2013).
  * (17) X.-J. Liu and H. Hu, Inhomogeneous Fulde-Ferrell superfluidity in spin-orbit-coupled atomic Fermi gases, Phys. Rev. A 87, 051608(R) (2013).
  * (18) L. He, H. Hu, and X.-J. Liu, Realizing Fulde-Ferrell Superfluids via a Dark-State Control of Feshbach Resonances, Phys. Rev. Lett. 120, 045302 (2018).
  * (19) I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * (20) S. Giorgini, L.P. Pitaevskii, S. Stringari, Theory of ultracold atomic Fermi gases. Rev. Mod. Phys. 80, 1215 (2008).
  * (21) G. Veeravalli, E. Kuhnle, P. Dyke, and C. J. Vale, Bragg Spectroscopy of a Strongly Interacting Fermi Gas, Phys. Rev. Lett. 101, 250403 (2008).
  * (22) S. Hoink, P. Dyke, M. G. Lingham, J. J. Kinnunen, G. M. Bruun and C. J. Vale, Goldstone mode and pair-breaking excitations in atomic Fermi superfluid, Nat. Phys. 13, 943 (2017).
  * (23) X. Li, X. Luo, S. Wang, K. Xie, X.-P. Liu, H. Hu, Y.-A. Chen, X.-C. Yao, J.-W. Pan, Second sound attenuation near quantum criticality, Science 375, 528 (2022).
  * (24) P. Dyke, S. Musolino, H. Kurkjian, D. J. M. Ahmed-Braun, A. Pennings, I. Herrera, S. Hoinka, S. J. J. M. F. Kokkelmans, V. E. Colussi, and C. J. Vale, Higgs oscillations in a unitary Fermi superfluid, arXiv:2310.03452v1.
  * (25) R. Senaratne, D. C.-Cavazos, S. Wang, F. He, Y.-T. Chang, A. Kafle, H. Pu, X.-W. Guan, R. G. Hulet, Spin-charge separation in a one-dimensional Fermi gas with tunable interactions, Science 376, 1305 (2022).
  * (26) P. W. Anderson, Random-phase approximation in the superconductivity, Phys. Rev. 112, 1900 (1958).
  * (27) X.-J. Liu, H. Hu, A. Minguzzi, and M. P. Tosi, Collective oscillations of a confined Bose gas at finite temperature in the random-phase approximation, Phys. Rev. A 69, 043605 (2004).
  * (28) K. V. Samokhin, Goldstone modes in Larkin-Ovchinnikov-Fulde-Ferrell superconductors, Phys. Rev. B 81, 224507 (2010).
  * (29) Z. Huang, C. S. Ting, J.-X. Zhu, and S.-Z. Lin, Gapless Higgs mode in the Fulde-Ferrell-Larkin-Ovchinnikov state of a superconductor, Phys. Rev. B 105, 014502 (2022).
  * (30) M. O. J. Heikkinen and P. Törmä, Collective modes and the speed of sound in the Fulde-Ferrell-Larkin-Ovchinnikov state, Phys. Rev. A 83, 053630 (2011).
  * (31) Z. Gao, L. He, H. Zhao, S.-G. Peng, and P. Zou, Dynamic structure factor of one-dimensional Fermi superfluid with spin-orbit coupling, Phys. Rev. A 107, 013304 (2023).
  * (32) H. Zhao, X. Yan, S.-G. Peng, and P. Zou, Dynamic structure factor of two-dimensional Fermi superfluid with Rashba spin-orbit coupling, Phys. Rev. A 108, 033309 (2023).


11institutetext: Nanjing University of Science and Technology, Nanjing,
210094, Jiangsu, China 11email<EMAIL_ADDRESS>
22institutetext: Nanjing University of Chinese Medicine, Nanjing, 210023,
Jiangsu, China

# TCMD: A Traditional Chinese Medicine QA Dataset for Evaluating Large
Language Models

Ping Yu 11 Kaitao Song 11 Fengchen He 11 Ming Chen 22 Jianfeng Lu 11

###### Abstract

The recently unprecedented advancements in Large Language Models (LLMs) have
propelled the medical community by establishing advanced medical-domain
models. However, due to the limited collection of medical datasets, there are
only a few comprehensive benchmarks available to gauge progress in this area.
In this paper, we introduce a new medical question-answering (QA) dataset that
contains massive manual instruction for solving Traditional Chinese Medicine
examination tasks, called TCMD. Specifically, our TCMD collects massive
questions across diverse domains with their annotated medical subjects and
thus supports us in comprehensively assessing the capability of LLMs in the
TCM domain. Extensive evaluation of various general LLMs and medical-domain-
specific LLMs is conducted. Moreover, we also analyze the robustness of
current LLMs in solving TCM QA tasks by introducing randomness. The
inconsistency of the experimental results also reveals the shortcomings of
current LLMs in solving QA tasks. We also expect that our dataset can further
facilitate the development of LLMs in the TCM area.

###### Keywords:

Traditional Chinese Medicine LLM Evaluation Multiple-Choice Question Answer

##  1 Introduction

In recent years, medical research has emerged as a frontier field at the
intersection of science, medicine, and technological innovation, employing
state-of-the-art technologies to extend our understanding of health and
disease. Recently, driven by the success of large language models (LLMs) [8,
22, 12, 23], building LLM-based medical applications [21, 27, 24] has drawn
enormous concerns. To evaluate these models and applications, various medical
datasets [11, 34, 26] are introduced, which promotes the prosperity and
development of AI in healthcare.

Traditional Chinese Medicine (TCM) remains a crucial and fascinating field in
Chinese medicine. So, to advance the integration of TCM with LLMs, researchers
have collected or generated extensive TCM corpora and trained LLMs
specifically for TCM [29, 36, 31, 24]. Additionally, several Chinese LLM
providers have incorporated TCM data into their training datasets [2, 12] to
ensure these models are well-aligned with TCM knowledge. However, few of the
current datasets are designed to evaluate the abilities of LLMs in TCM as some
of them do not include samples from TCM, while some treat samples from TCM as
part of their comprehensive evaluation. Another important point is that most
LLMs or LLM-powered applications in the TCM domain rely on subjective criteria
such as usability and smoothness to assess their alignment with human
expectations. However, these subjective metrics are labour-intensive and often
lack consistency. Therefore, there is a pressing need for a new dataset that
enables objective evaluation to address this gap.

As [5] suggested, it is a natural option to use problems from standard
examinations to objectively assess medical LLMs. For example, [15, 14, 19]
propose to use questions from the United States Medical Licensing Examination
(USMLE) to test the performance of LLMs. Some of the tested LLMs perform well
on these datasets, and a few of them even surpass the passing line for humans.

Following this principle, we propose building a multiple-choice QA (MCQA)
dataset called TCMD to evaluate LLMs’ performance in TCM objectively. Inspired
by [17]and [26], we collect multiple-choice problems related to the Chinese
National Medical Licensing Examination For Traditional Chinese Medicine and
their corresponding explanations. Subsequently, we filter and organize the
questions according to the official examination manual to ensure comprehensive
coverage of all subjects mentioned in the instructions. This process
guarantees that the distribution of questions across various medical topics
aligns with the guidelines. Furthermore, we introduce two types of multiple-
choice problems not included in the two benchmarks. With this new dataset, we
benchmark the performance of general LLMs, common medical LLMs, and TCM-
domain-specific LLMs. We also propose to evaluate the robustness of LLMs by
introducing randomness in the test, i.e., testing LLMs with the question whose
options are shuffled. Surprisingly, the consistency of the predicted answers
is quite poor, even for LLMs with high scores in normal evaluation. More
details of the result will be discussed in Section 4.

Overall, the main contributions and findings of this paper can be summarized
as:

  1. 1.

We build a new dataset to evaluate the capability of different LLMs in solving
TCM tasks. It contains 4 types of multiple-choice questions and their
explanations in the style of the Chinese National Medical Licensing
Examination. The questions are filtered and organized under the guidance of
the official manual instruction.

  2. 2.

We conducted extensive analysis to investigate different LLMs and found that
the general LLMs perform better than medical and TCM LLMs on average.
Moreover, the consistency of answers predicted by LLMs when facing questions
with shuffled options is far from satisfying.

  3. 3.

Moreover, we investigate the ensemble of predictions for questions with
shuffled options using a voting mechanism and discover that this ensemble
strategy can enhance final performance under certain conditions.

##  2 Related Work

###  2.1 LLMs for Medical Application

In recent years, there has been an increasing interest in the application of
LLMs in the field of medicine. Researchers have explored various approaches to
leverage LLMs for medical tasks such as biomedical information extraction [13,
25, 32], clinical text analysis [32, 6], medical queries [33, 10], and so on.
Certain researchers take an additional stride by directly replacing the doctor
with LLMs [30, 16]. Particularly, in a double-blind trial, MedGPT, proposed by
Medlinker, demonstrated a medical concordance rate of 96% with attending
physicians from Three-A hospital111https://www.medlinker.com/news/201. Also,
Med-PALM 2 [21] becomes the first LLM to pass the United States Medical
Licensing Examination style questions in MultiMedQA [20] with an expert-level
score. Through such a plethora of studies, we can envision that in the near
future, LLMs will be a great help to health workers.

###  2.2 Medical Evaluation Datasets for LLMs

Medical evaluation datasets for LLM are hard to construct as they require
well-designed samples with diverse formats and rich medical knowledge. He [11]
proposes to collect medical questions from online health consultancy websites
to construct medical QA datasets. The dataset is unsatisfying in terms of
professionalism and fluency because the conversations on online platforms are
usually not conventional. To overcome the problem, Li [15] proposes the first
large-scale open-source Chinese medical MCQA dataset, i.e. MLEC-QA, sourced
from the Chinese National Medical Licensing Examination. Liu [17] collects
problems from the same source and performs a more detailed categorization,
making it possible to measure in multiple dimensions. To better leverage the
existing datasets, [20] collects several medical datasets covering medical
questions from multiple countries to assess LLMs in medical question
answering.

Most of these datasets collect as many problems as possible to evaluate
medical abilities comprehensively. However, knowledge of different subjects
holds varying significance in real-world applications. To address it, our
dataset contains questions from every required medical subject in the national
exam, and we also balance the distribution of questions of different subjects
according to their official manual.

##  3 Dataset Construction

###  3.1 Manual Instruction of CNMLE

The manual instruction of CNMLE gives detailed information about the CNMLE,
including question types, exam content, exam schedule, etc. It is curated by
medical and educational experts to evaluate students’ mastery of essential
knowledge required for real-world medical practice. As a result, it can be the
guiding principle to construct an evaluation dataset.

##### Question Types

According to the part for TCM, there are 4 types of questions:

  * •

A1\A2: a single multiple-choice problem with a statement or a clinical case.

  * •

A3: several questions with a shared stem consisting of a clinical history.

  * •

B1: several questions with shared options containing interrelated concepts.

For convenience, we set the number of sub-problems for A3 and B1 problems as 3
and 2. English examples for different types can be found in the dataset.

Figure 1: The detailed exam content of the CNMLE for TCM. The mark * indicates
that the subject is not required for those with apprenticeships or
demonstrable expertise.

##### Exam Content

The instruction indicates that 5 medical realms and 18 subjects should be
involved. A detailed list of them is shown in Figure 1, and so is the
percentage of questions related to each realm. It’s noticeable that the
Clinical Practice of TCM accounts for half of the questions, which shows its
importance in real-world diagnosis. The exam contains 600 questions in total.

###  3.2 Data Collecting and Processing

##### Problems Collection

Most of the questions are from two sources: reference books and TCM
enthusiasts. Many of them are from actual examinations in different years. We
also collect an explanation for each problem to consult the reasoning of LLMs.

##### Manual Verification and Processing

Most of these materials are in PDF format, so we first use optical character
recognition to convert them to text files. Then, the files are verified by
humans and programs to ensure their correctness. After verification, we de-
duplicate the questions and group them based on which exam unit they belong
to. Finally, medical experts are invited to check the examples.

##### Balance of Dataset

To ensure that the dataset covers all the subjects and the test set follows
the distribution from the manual instruction, we first invite some experts to
label the subject of each collected problem. We drop or extend problems
according to the gap between the real number and the target number of problems
for different realms.

##### Data Statistics

Finally, we collected 2851 questions for the train set and 600 questions for
the test set. The number of questions of different types for different
subjects is shown in Table 1.

Table 1: The number of questions of different types for each subject grouped by realm. The first number indicates the train set, and the second is for the test set.  Subject | A1 | A2 | A3 | B1 | Total
---|---|---|---|---|---
Fundamental Theory of TCM | 259/25 | 3/0 | 0/0 | 87/10 | 349/35
Diagnostic Theory of TCM | 197/17 | 136/6 | 45/0 | 131/12 | 509/35
Chinese Herbal Medicine | 150/15 | 13/6 | 0/0 | 85/14 | 248/35
Prescription Formulation Studies for TCM | 250/19 | 126/8 | 15/0 | 141/8 | 532/35
Huangdi’s Internal Classic | 4/5 | 0/0 | 0/0 | 0/0 | 4/5
Treatise on Febrile Diseases | 5/1 | 0/4 | 0/0 | 0/0 | 5/5
Synopsis of the Golden Chamber | 4/1 | 0/4 | 0/0 | 0/0 | 4/5
epidemic febrile diseases | 0/2 | 0/3 | 0/0 | 1/0 | 1/5
Internal Medicine of TCM | 4/18 | 0/27 | 0/27 | 8/18 | 12/90
Surgical Medicine of TCM | 59/11 | 16/15 | 0/12 | 8/12 | 83/50
Gynecology of TCM | 74/10 | 47/15 | 0/15 | 11/10 | 132/50
Pediatrics of TCM | 54/10 | 41/15 | 0/15 | 3/10 | 98/50
Acupuncture and Moxibustion | 114/12 | 60/18 | 0/18 | 37/12 | 211/60
Internal Medicine | 55/10 | 21/20 | 0/15 | 36/14 | 112/59
Fundamentals of Diagnostic Medicine | 209/32 | 69/1 | 0/0 | 76/8 | 354/41
Infectious Disease Epidemiology | 62/12 | 3/0 | 0/0 | 19/8 | 84/20
Medical Ethics | 35/6 | 0/0 | 0/0 | 19/4 | 54/10
Health Regulations and Legislation | 45/6 | 1/0 | 0/0 | 13/4 | 59/10
Table 2: The performance of each model on the test set. The FoT, CoT, CPoT, IWM, MH in the header are short for Foundation of TCM, Classics of TCM, Clinical Practice of TCM, Integrated Western Medicine, Medical Humanities respectively. In the method column, "ICL" stands for In-Context Learning, and "COT" represents Chain-of-Thought. Models marked with an asterisk (*) are closed-source commercial models. A dash (-) in the score indicates that the model did not generate answers in the required format. The $\uparrow$ signifies the realm in which the model excels most out of the five domains, whereas $\downarrow$ indicates its weakest performance. Model Name | Domain | Method | REALM | Total
---|---|---|---|---
FoT | CoT | CPoT | IWM | MH
ChatGPT4* | General | ICL | 0.550 | 0.350 | 0.563 | 0.667 | 0.850 | 0.583
ICL + COT | 0.543 | 0.350 | 0.593 | 0.550 | 0.650 | 0.567
GPT-35-turbo-instruct* | General | ICL | 0.357 | 0.250 | 0.410 | 0.475 | 0.750 | 0.417
ICL + COT | 0.350 | 0.000 | 0.310 | 0.250 | 0.550 | 0.305
Moonshot-v1-8k* | General | ICL | 0.693 | 0.500 | 0.713 | 0.708 | 0.600 | 0.697
ICL + COT | 0.793 | 0.550 | 0.780 | 0.758 | 0.700 | 0.768
AquilaChat2-34B | General | ICL | 0.421 | 0.300 | 0.407 | 0.483 | 0.650 | 0.430
ICL + COT | 0.271 | 0.150 | 0.310 | 0.317 | 0.200 | 0.293
AquilaChat2-7B | General | ICL | 0.293 | 0.350 | 0.287 | 0.375 | 0.350 | 0.310
ICL + COT | 0.236 | 0.200 | 0.263 | 0.217 | 0.250 | 0.245
Baichuan2-13B-Base | General | ICL | 0.657 | 0.600 | 0.607 | 0.592 | 0.650 | 0.617
ICL + COT | 0.657 | 0.500 | 0.617 | 0.642 | 0.550 | 0.625
Baichuan2-7B-Base | General | ICL | 0.536 | 0.400 | 0.547 | 0.617 | 0.350 | 0.547
ICL + COT | 0.600 | 0.550 | 0.520 | 0.617 | 0.550 | 0.560
Baichuan-13B-Base | General | ICL | 0.529 | 0.300 | 0.483 | 0.492 | 0.550 | 0.492
ICL + COT | 0.486 | 0.450 | 0.410 | 0.458 | 0.750 | 0.450
Baichuan-7B | General | ICL | 0.371 | 0.250 | 0.380 | 0.375 | 0.350 | 0.372
ICL + COT | 0.314 | 0.450 | 0.317 | 0.442 | 0.400 | 0.348
ChatGLM3-6B | General | ICL | 0.350 | 0.150 | 0.407 | 0.383 | 0.300 | 0.377
ICL + COT | 0.464 | 0.250 | 0.420 | 0.400 | 0.300 | 0.417
ChatGLM2-6B | General | ICL | 0.421 | 0.200 | 0.403 | 0.475 | 0.400 | 0.415
ICL + COT | 0.450 | 0.200 | 0.480 | 0.542 | 0.550 | 0.478
ChatGLM-6B | General | ICL | 0.321 | 0.250 | 0.250 | 0.350 | 0.400 | 0.292
ICL + COT | 0.257 | 0.350 | 0.267 | 0.275 | 0.250 | 0.268
InternLM2-Chat-20B | General | ICL | 0.614 | 0.400 | 0.653 | 0.650 | 0.600 | 0.633
ICL + COT | 0.671 | 0.500 | 0.657 | 0.650 | 0.700 | 0.655
InternLM2-Chat-7B | General | ICL | 0.557 | 0.300 | 0.660 | 0.500 | 0.450 | 0.585
ICL + COT | 0.550 | 0.250 | 0.593 | 0.525 | 0.650 | 0.560
Qwen-14B-Chat | General | ICL | 0.629 | 0.150 | 0.570 | 0.667 | 0.700 | 0.593
ICL + COT | 0.629 | 0.350 | 0.563 | 0.675 | 0.700 | 0.598
Qwen-7B-Chat | General | ICL | 0.400 | 0.200 | 0.443 | 0.467 | 0.450 | 0.430
ICL + COT | 0.464 | 0.500 | 0.467 | 0.467 | 0.650 | 0.473
Vicuna-7B | General | ICL | 0.264 | 0.050 | 0.227 | 0.292 | 0.200 | 0.242
ICL + COT | 0.157 | 0.100 | 0.163 | 0.142 | 0.200 | 0.157
ChatMed-Consult | Medical | ICL | 0.157 | 0.200 | 0.177 | 0.133 | 0.250 | 0.167
ICL + COT | 0.050 | 0.000 | 0.030 | 0.008 | 0.150 | 0.033
ShenNong-TCM-LLM | TCM | ICL | 0.129 | 0.300 | 0.100 | 0.100 | 0.300 | 0.120
ICL + COT | 0.064 | 0.050 | 0.030 | 0.050 | 0.100 | 0.045
Huatuo-Llama-Med-Chinese | TCM | ICL | 0.229 | 0.350 | 0.127 | 0.117 | 0.250 | 0.160
ICL + COT | - | - | - | - | - | -

##  4 Experiments and Analysis

To make a thorough evaluation, we test various models on our datasets. The
details of the experiment will be shown in this section.

###  4.1 Experimental Setup

##### Baseline

We choose models of three different categories.

  * •

General LLMs: ChatGPT(ChatGPT4, GPT-35-turbo-instruct) [18],
Kimi(Moonshot-v1-8k) 222https://platform.moonshot.cn/, ChatGLM [8], Baichuan
[2], QwenChat [22], InternLMChat [12], AquilaChat [1] and Vicuna [35].

  * •

Common Medical LLMs: ChatMed [36]

  * •

LLMs for TCM: ShenNong-TCM-LLM333https://github.com/michael-wzhu/ShenNong-TCM-
LLM, Huatuo-Llama-Med-Chinese [24]

Most models are open-source and can be found on the provided URL or
Huggingface 444https://huggingface.co/.

##### Prompt, Methods and Settings

We hand-craft some prompts, while the others are from the model providers’
evaluation code. These prompts can be found in the dataset. For prompting
methods, we employ In-Context-Learning [4] and Chain-of-Thought [28].
Specifically, we retrieve 3 examples for a demonstration from the train set
using the Rank-BM25 555https://github.com/dorianbrown/rank_bm25. We also think
splitting A3 and B1 questions into single sub-questions is against their
purpose, so we feed the LLMs with the whole question. As for the generation
parameters, we use the default ones except that we set do_sampe to False to
avoid irreproducible results. We set max_new_tokens to 128. For the same
purpose, we set seed to 0 for ChatGPT and Kimi.

###  4.2 Main Result

We use accuracy as our evaluation metric. The detailed scores for each model
are listed in Table 2. Despite the fact that the medical and TCM LLMs used
medical or TCM corpus during training, the overall performances of general
LLMs are much better. We think the tuning process of medical LLMs and TCM LLMs
only involves textual QA tasks, which may result in catastrophic forgetting
[9] of their abilities to answer multiple-choice questions. Among these LLMs,
Moonshot-v1-8k using chain-of-thought earns the highest score. Among those
open-source models, InternLM2-Chat-20B using chain-of-thought archives the
best performance. Theoretically, if the testee answers correctly to over 60%
of the questions, he will pass the examination. So, Moonshot-v1-8k,
Baichuan2-13b-Base and InternLM2-Chat-20B are the three passers of our
evaluation. The training data of Baichuan2 and InternLM2 includes medical
literature according to their technical report, and we think it contributes to
their outstanding performance. Another noteworthy point is that ChatGPT4, one
of the most powerful LLMs in the world, does not achieve a passing score, and
the reason might be the lack of TCM knowledge in the training phase.

###  4.3 Analysis

We will analyze our results from the following four perspectives.

##### The performance in each realm

Among all the 5 realms, most models perform terribly in Classics of TCM. We
think the reason for the poor performance is that classical Chinese literature
is not common in the training corpus but is ubiquitous in this realm. Although
Classics of TCM only account for 3.33% of the questions, it plays an important
role in learning TCM, and the Chinese people have a strong passion for
classical Chinese.

##### The performance on different question types

Statistics about performance on different types are shown in Table 3. To avoid
redundancy, those results whose total accuracy is no less than 0.5 are
analyzed. The statistics show that most of the LLMs perform poorly on type B1
and type A3. Type B1 questions usually contain confusing concepts, and some
medical experts may even make mistakes about these concepts. Type A3 questions
contain complex reasoning processes. So, the abilities of LLMs for concept
discrimination and reasoning in TCM should be strengthened.

##### The birth time of models

For each series of models, like ChatGLM and Baichuan, we can see that the
newer versions of models are generally better than the old ones. Particularly,
Baichuan2-7B-Base even beats Baichuan-13B-Base despite the huge gap in
parameter size. The phenomenon proves that medical applications hold a
significant strategic position among various LLM providers.

##### Is Chain-of-Thought a universal remedy?

Chain-of-Thought (COT) is a powerful and well-validated reasoning method that
can enhance the logical reasoning capabilities of LLMs. However, our
experimental results show that using COT did not improve accuracy for several
models. We believe this is because most models have not been trained on data
from the TCM domain, leading to a lack of understanding of many TCM concepts.
As a result, using COT may conflict with the models’ internal knowledge,
reducing their performance. Therefore, COT cannot be considered a universal
solution for enhancing model performance in the TCM domain.

Table 3: The performance of each model on different types of questions.  Model | Method | A1 | A2 | A3 | B1
---|---|---|---|---|---
ChatGPT4* | ICL | 0.580 | 0.613 | 0.569 | 0.569
ICL + COT | 0.561 | 0.465 | 0.647 | 0.618
moonshot-v1-8k* | ICL | 0.741 | 0.746 | 0.686 | 0.590
ICL + COT | 0.807 | 0.768 | 0.686 | 0.771
Baichuan2-13B-Base | ICL | 0.651 | 0.725 | 0.588 | 0.479
ICL + COT | 0.656 | 0.669 | 0.539 | 0.597
Baichuan2-7B-Base | ICL | 0.599 | 0.662 | 0.520 | 0.375
ICL + COT | 0.594 | 0.641 | 0.510 | 0.465
InternLM2-Chat-20B | ICL | 0.675 | 0.683 | 0.608 | 0.542
ICL + COT | 0.646 | 0.662 | 0.627 | 0.681
InternLM2-Chat-7B | ICL | 0.608 | 0.634 | 0.686 | 0.431
ICL + COT | 0.580 | 0.585 | 0.569 | 0.500
Qwen-14B-Chat | ICL | 0.689 | 0.599 | 0.588 | 0.451
ICL + COT | 0.679 | 0.683 | 0.412 | 0.528
Table 4: The statistics of accuracy and consistency for predictions to shuffled variants of the test set. Model | Accuracy | Consistency
---|---|---
| Mean(STD) | Min | Max | Mean(STD) | Min | Max
Baichuan2-13B-Base | 0.593(0.016) | 0.560 | 0.633 | 0.767(0.224) | 0.200 | 1.000
Baichuan2-7B-Base | 0.540(0.015) | 0.505 | 0.580 | 0.699(0.221) | 0.200 | 1.000
Baichuan-13B-Base | 0.481(0.018) | 0.442 | 0.530 | 0.648(0.212) | 0.210 | 1.000
ChatGLM3-6B | 0.385(0.018) | 0.335 | 0.437 | 0.604(0.225) | 0.200 | 1.000
ChatGLM2-6B | 0.397(0.019) | 0.360 | 0.438 | 0.584(0.259) | 0.200 | 1.000
Qwen-14B-Chat | 0.603(0.017) | 0.570 | 0.645 | 0.756(0.223) | 0.200 | 1.000
Qwen-7B-Chat | 0.444(0.015) | 0.400 | 0.478 | 0.623(0.231) | 0.200 | 1.000
InternLM2-Chat-7B | 0.592(0.012) | 0.570 | 0.628 | 0.795(0.220) | 0.200 | 1.000
Figure 2: The accuracy of LLMs under different thresholds of consistency. The
dashed line is the mean accuracy of the 120 variants from Table 4.

###  4.4 Robustness Evaluation

Some recent studies [3] also indicate the LLM could be sensitive for the
training order (e.g., a model trained on “A is B” cannot be generalized to the
reverse pattern “B is A”). Therefore, it is necessary to validate whether LLMs
are robust in real-world applications. To this end, we also conduct
experiments for robustness evaluation by using a simple strategy.

####  4.4.1 Test Strategy

To validate the robustness of LLMs, we review our preparation process for an
exam during which we may exercise countless questions. To ensure that we grasp
the question instead of simply memorizing it, we may shuffle the options and
check the consistency of our answers. We consider this strategy to be useful
for evaluating LLMs. So, we propose to shuffle the options and check whether
the LLM gives a consistently good performance. To begin with, we define the
consistency for question $q$ as

$\operatorname{cons}(q)=\frac{1}{n}\operatorname{max}(\operatorname{freq}(ans,pred),ans=A,B,C,D,E),$
(1)

where $n$ is the number of variants for the question, $pred$ is the $n$
predictions made by the LLM and $\operatorname{freq}(ans,pred)$ gives the
frequency of $ans$ among $pred$.

To minimize the influence of randomness, we collect the results of all
possible mappers, i.e., the full permutations of the options. As we have 5
options for each question, we will have 120 shuffled variants for 1 question.
Limited by resource, we only test 8 open-source well-performed models whose
size is below 14B, which can be found in 4 and the prompt method is ICL.

####  4.4.2 Evaluation Result and Analysis

We present and analyze the test results from the perspectives of accuracy and
consistency.

##### Accuracy

The statistical figures of models’ accuracy across all shuffled variants of
the test set are shown in Tabel 4. The gap between mean accuracy and the
accuracy in Table 2 is narrow, and the standard diversion of accuracy remains
small. These figures indicate that the overall performance of LLMs is stable.

##### Consistency

We calculate the statistics of the consistency across the test set, which is
shown in Table 4. As the figures indicate, even when facing two nearly
identical questions, LLMs still have an unignorable chance to output different
answers.

####  4.4.3 Analysis on Ensemble Methods

A common and robust way to deal with predictions with uncertainty is to
integrate multiple predictions, i.e., ensemble [7]. We further investigate how
the inconsistency of LLMs may affect this strategy by integrating the results.
We take the most frequently occurring result among the 120 predicted answers
as the final answer if the consistency passes the threshold. Otherwise, we
consider the LLMs failed to generate the answer. The relationship between the
threshold and the accuracy of the integrated result is shown in Figure 2. For
all models, the accuracy drops as the threshold rises. Particularly, the
performance becomes unacceptable when the threshold goes high, which makes
their robustness questionable. Interestingly, the integrated accuracy is
higher than the mean accuracy from Table 4, which suggests that the
integration may boost the performance if the threshold is very low.

####  4.4.4 Summary

Although the overall performance of most LLMs remains stable when tested with
shuffled variants of questions, the consistency of all the LLMs is not so
promising. Fortunately, our experiments suggest that ensemble methods can
alleviate this problem and may even boost the performance of LLMs. However, we
should still think twice about adopting LLMs when faced with risk-averse tasks
in complex and diverse environments, like surgery and financial systems.

##  5 Conclusion

In this paper, we propose a new dataset to systemically evaluate LLMs for
their ability in TCM with MCQA questions under the guidance of the official
examination manual instruction for CNMLE. The dataset contains questions from
various realms and medical subjects. We also propose to evaluate the
robustness of LLMs by the consistency of predictions when facing questions
whose options are shuffled. We hope the dataset and our evaluation method will
foster the advancement of LLMs in the field of TCM.

## References

  * [1] BAAI: Aquilachat. https://huggingface.co/BAAI (2023)
  * [2] Baichuan: Baichuan 2: Open large-scale language models. arXiv preprint arXiv:2309.10305 (2023), https://arxiv.org/abs/2309.10305
  * [3] Berglund, L., Tong, M., Kaufmann, M., Balesni, M., Stickland, A.C., Korbak, T., Evans, O.: The Reversal Curse: LLMs trained on "A is B" fail to learn "B is A" (Sep 2023), http://arxiv.org/abs/2309.12288, arXiv:2309.12288 [cs]
  * [4] Brown, T., Mann, B., Ryder, N., Subbiah, M., Kaplan, J.D., Dhariwal, P., Neelakantan, A., Shyam, P., Sastry, G., Askell, A., et al.: Language models are few-shot learners. Advances in neural information processing systems 33, 1877–1901 (2020)
  * [5] Chang, Y., Wang, X., Wang, J., Wu, Y., Zhu, K., Chen, H., Yang, L., Yi, X., Wang, C., Wang, Y., et al.: A survey on evaluation of large language models. arXiv preprint arXiv:2307.03109 (2023)
  * [6] Chervenak, J., Lieman, H., Blanco-Breindel, M., Jindal, S.: The promise and peril of using a large language model to obtain clinical information: Chatgpt performs strongly as a fertility counseling tool with limitations. Fertility and Sterility (2023)
  * [7] Dietterich, T.G.: Ensemble methods in machine learning. In: International workshop on multiple classifier systems. pp. 1–15. Springer (2000)
  * [8] Du, Z., Qian, Y., Liu, X., Ding, M., Qiu, J., Yang, Z., Tang, J.: Glm: General language model pretraining with autoregressive blank infilling. In: Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers). pp. 320–335 (2022)
  * [9] French, R.M.: Catastrophic forgetting in connectionist networks. Trends in cognitive sciences 3(4), 128–135 (1999)
  * [10] Hashemi, H., Zhuang, Y., Kothur, S.S.R., Prasad, S., Meij, E., Croft, W.B.: Dense retrieval adaptation using target domain description. In: Proceedings of the 2023 ACM SIGIR International Conference on Theory of Information Retrieval. pp. 95–104 (2023)
  * [11] He, J., Fu, M., Tu, M.: Applying deep matching networks to chinese medical question answering: a study and a dataset. BMC medical informatics and decision making 19(2), 91–100 (2019)
  * [12] InternLMTeam: Internlm: A multilingual language model with progressively enhanced capabilities. https://github.com/InternLM/InternLM (2023)
  * [13] Kartchner, D., Ramalingam, S., Al-Hussaini, I., Kronick, O., Mitchell, C.: Zero-shot information extraction for clinical meta-analysis using large language models. In: The 22nd Workshop on Biomedical Natural Language Processing and BioNLP Shared Tasks. pp. 396–405 (2023)
  * [14] Kung, T.H., Cheatham, M., Medenilla, A., Sillos, C., De Leon, L., Elepaño, C., Madriaga, M., Aggabao, R., Diaz-Candido, G., Maningo, J., et al.: Performance of chatgpt on usmle: Potential for ai-assisted medical education using large language models. PLoS digital health 2(2), e0000198 (2023)
  * [15] Li, J., Zhong, S., Chen, K.: Mlec-qa: A chinese multi-choice biomedical question answering dataset. In: Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing. pp. 8862–8874 (2021)
  * [16] Li, Y., Li, Z., Zhang, K., Dan, R., Jiang, S., Zhang, Y.: Chatdoctor: A medical chat model fine-tuned on a large language model meta-ai (llama) using medical domain knowledge. Cureus 15(6) (2023)
  * [17] Liu, J., Zhou, P., Hua, Y., Chong, D., Tian, Z., Liu, A., Wang, H., You, C., Guo, Z., Zhu, L., et al.: Benchmarking large language models on cmexam–a comprehensive chinese medical exam dataset. arXiv preprint arXiv:2306.03030 (2023)
  * [18] OpenAI: Chatgpt: Language model. https://www.openai.com/chatgpt (2023)
  * [19] Sharma, P., Thapa, K., Dhakal, P., Upadhaya, M.D., Adhikari, S., Khanal, S.R.: Performance of chatgpt on usmle: Unlocking the potential of large language models for ai-assisted medical education. arXiv preprint arXiv:2307.00112 (2023)
  * [20] Singhal, K., Azizi, S., Tu, T., Mahdavi, S.S., Wei, J., Chung, H.W., Scales, N., Tanwani, A., Cole-Lewis, H., Pfohl, S., et al.: Large language models encode clinical knowledge. arXiv preprint arXiv:2212.13138 (2022)
  * [21] Singhal, K., Tu, T., Gottweis, J., Sayres, R., Wulczyn, E., Hou, L., Clark, K., Pfohl, S., Cole-Lewis, H., Neal, D., et al.: Towards expert-level medical question answering with large language models. arXiv preprint arXiv:2305.09617 (2023)
  * [22] Team, Q.: Qwen technical report. arXiv preprint arXiv:2309.16609 (2023)
  * [23] Touvron, H., Martin, L., Stone, K., Albert, P., Almahairi, A., Babaei, Y., Bashlykov, N., Batra, S., Bhargava, P., Bhosale, S., et al.: Llama 2: Open foundation and fine-tuned chat models. arXiv preprint arXiv:2307.09288 (2023)
  * [24] Wang, H., Liu, C., Xi, N., Qiang, Z., Zhao, S., Qin, B., Liu, T.: Huatuo: Tuning llama model with chinese medical knowledge (2023)
  * [25] Wang, Q., Gao, Z., Xu, R.: Exploring the in-context learning ability of large language model for biomedical concept linking. arXiv preprint arXiv:2307.01137 (2023)
  * [26] Wang, X., Chen, G.H., Song, D., Zhang, Z., Chen, Z., Xiao, Q., Jiang, F., Li, J., Wan, X., Wang, B., et al.: Cmb: A comprehensive medical benchmark in chinese. arXiv preprint arXiv:2308.08833 (2023)
  * [27] Wang, Z., Li, R., Dong, B., Wang, J., Li, X., Liu, N., Mao, C., Zhang, W., Dong, L., Gao, J., et al.: Can llms like gpt-4 outperform traditional ai tools in dementia diagnosis? maybe, but not today. arXiv preprint arXiv:2306.01499 (2023)
  * [28] Wei, J., Wang, X., Schuurmans, D., Bosma, M., Xia, F., Chi, E., Le, Q.V., Zhou, D., et al.: Chain-of-thought prompting elicits reasoning in large language models. Advances in Neural Information Processing Systems 35, 24824–24837 (2022)
  * [29] Wei Zhu, W.Y., Wang, X.: Shennong-tcm: A traditional chinese medicine large language model. https://github.com/michael-wzhu/ShenNong-TCM-LLM (2023)
  * [30] Xiong, H., Wang, S., Zhu, Y., Zhao, Z., Liu, Y., Wang, Q., Shen, D.: Doctorglm: Fine-tuning your chinese doctor is not a herculean task. arXiv preprint arXiv:2304.01097 (2023)
  * [31] Yang, S., Zhao, H., Zhu, S., Zhou, G., Xu, H., Jia, Y., Zan, H.: Zhongjing: Enhancing the chinese medical capabilities of large language model through expert feedback and real-world multi-turn dialogue. arXiv preprint arXiv:2308.03549 (2023)
  * [32] Yang, X., Chen, A., PourNejatian, N., Shin, H.C., Smith, K.E., Parisien, C., Compas, C., Martin, C., Costa, A.B., Flores, M.G., et al.: A large language model for electronic health records. NPJ Digital Medicine 5(1), 194 (2022)
  * [33] Zakka, C., Chaurasia, A., Shad, R., Hiesinger, W.: Almanac: Knowledge-grounded language models for clinical medicine. arXiv preprint arXiv:2303.01229 (2023)
  * [34] Zhang, N., Chen, M., Bi, Z., Liang, X., Li, L., Shang, X., Yin, K., Tan, C., Xu, J., Huang, F., et al.: Cblue: A chinese biomedical language understanding evaluation benchmark. arXiv preprint arXiv:2106.08087 (2021)
  * [35] Zheng, L., Chiang, W.L., Sheng, Y., Zhuang, S., Wu, Z., Zhuang, Y., Lin, Z., Li, Z., Li, D., Xing, E., et al.: Judging llm-as-a-judge with mt-bench and chatbot arena. arXiv preprint arXiv:2306.05685 (2023)
  * [36] Zhu, W., Wang, X.: Chatmed: A chinese medical large language model. https://github.com/michael-wzhu/ChatMed (2023)