Chemical potential (in)dependence of hadron scatterings in the hadronic phase of QCD-like theories and its applications

We formulate a method for calculating the hadron-hadron scattering amplitudes at nonzero chemical potential ($\mu$) in the hadronic phase at zero temperature, where the baryon number symmetry remains to be violated. Although it is widely believed that the physical quantities do not change even if we turn on a small $\mu$ at zero temperature, the shape of correlation functions for a single hadron depends on $\mu$. Then, the dispersion relation of the single hadron is modified to $E({\bf p},\mu) = \sqrt{{\bf p}^2+m^2}-\mu n_{O}$. Here, $m$ and $n_O$ denote the hadron mass at $\mu=0$ and the quantum number, respectively. From this relation, it is possible that the effective mass of the hadron depends on $\mu$. We extend the HAL QCD method at $\mu=0$ to the case of $\mu \ne 0$, which allows us to extract the scattering phase shifts via the interaction potential. We have found that the interaction potential can depend on $\mu$ only through the effective mass while the scattering phase shifts, obtained by solving the Schr\"{o}dinger equation with the interaction potential, are independent of $\mu$. We also numerically analyze the S-wave scatterings of two pions with isospin $I=2$ and two scalar diquarks within the framework of QC$_{2}$D at nonzero quark chemical potential. While the lattice is not exactly set to zero temperature, the $\mu$-independence can be observed. Furthermore, we improve the results for the S-wave scatterings of two hadrons obtained above by taking the $\mu$-independence for granted. Thanks to the asymmetric property of the correlation functions for diquarks at $\mu\neq0$, we can access a long-$\tau$ regime and can reduce the systematic error coming from inelastic contributions.


Introduction
The nature of Quantum Chromo Dynamics (QCD) and hadrons at finite density and low temperature (T ) is still poorly understood despite the presence of various theoretical and experimental studies.It is widely believed that a zero-temperature phase transition occurs from the hadronic to the nuclear superfluid phase at a critical chemical potential (µ c ) where the baryon number density n B jumps from zero to the nuclear saturation density n 0 (see Figure 1, where µ q and T c denote the quark chemical potential and the pseudocritical temperature at µ q = 0, respectively.)Once n B becomes nonzero, therefore, the properties of the quantum vacuum |0⟩ will drastically change, which has relevance to the internal structure of neutron stars.In contrast to such a high-density regime, the physical observables in equilibrium must be unchanged below µ c since the chemical potential (µ) acts as a source of particle production and hence µ c corresponds to the threshold value of dynamical baryon production.Therefore, the quark (or baryon) number must be conserved below µ c .This salient property is called the Sliver Blaze phenomenon [1].
Turning to the first-principles calculations for QCD, the nonzero chemical potential regime is still highly challenging due to the difficulty of the so-called sign problem [2].In the conventional Monte Carlo method, which is based on importance sampling, the probability weight given by the Euclidean action becomes complex if we add the number operator to the QCD action.In particular, the low-temperature and high-density regime is extremely difficult to access by the simulation although several new approaches are known [3][4][5][6][7][8][9][10][11].According to an analytical study based on the lattice gauge theory [12][13][14], the µ-independence of the quark number density for 0 < µ q < m π /2 has been proven at the T → 0 limit, while it is not exactly proven but ⟨n q ⟩ = 0 is suggested for m π /2 < µ q < m N /3 after taking the ensemble average [15].On the other hand, the two-point function for hadrons, which is not a local quantity, has a µ-dependence if we calculate it using the QCD action with a small chemical potential; the chemical potential term breaks imaginary-time reversibility.It is analytically known how these correlation functions and the corresponding hadron spectra depend on µ.
This paper discusses the µ-independence of the hadron scattering parameters, which can be numerically calculated by multi-point correlation functions of hadrons.We reformulate a HAL QCD method at µ = 0 [16][17][18] to µ ̸ = 0 for studying hadron scatterings with sufficiently small chemical potential at T = 0 and show what underlies the µ-independence and how it appears.In this process, we will clarify which parts of our formulation for µ q < µ c become invalid if we would consider the superfluid phase ranging µ q > µ c .The HAL QCD method helps to extract the information on the S-matrix via hadron interaction potentials and has been successful at µ = 0.For instance, it provides the nuclear forces at large quark masses, which are qualitatively consistent with the phenomenological one such as the AV18 potential [16][17][18].It has also been applied to the studies of tetraquarks [19][20][21][22], pentaquarks [23][24][25], and dibaryons [26][27][28][29][30]. Furthermore, applications of this method to hadron resonances have been carried out recently [31][32][33][34].
In the process of formulating the HAL QCD method at µ ̸ = 0, we do not specify the number of colours in the SU(N c ) gauge theory.Thus, the formulation can be applied to 3-colour QCD, if we could overcome the sign problem and could generate gauge configurations of dense 3-colour QCD.At present, however, it is an extremely hard task.In the latter part of this paper, we examine the hadron scatterings in 2-colour QCD (QC 2 D) at finite µ q to demonstrate our predictions of the µ-dependence.Dense QC 2 D coupled with an even number of quarks is free from the sign problem even for the conventional Monte Carlo simulations.The studies in such theories also serve as clues to understanding QCD in a medium.As in dense 3-colour QCD, it is expected that a phase transition from the hadronic to the superfluid phase occurs in dense 2-colour QCD at sufficiently low temperatures.In the case of QC 2 D, the baryon consists of two quarks, so that it is a bosonic particle.The 2-flavour QC 2 D model has an extended flavour symmetry, so-called Pauli-Gürsey symmetry [35,36], then the nucleon (diquark) mass m N becomes equivalent to the pion mass (m π ) at µ q = 0. Therefore, the critical chemical potential µ c is expected as µ c = m π /2.Actually, the emergence of the superfluid phase in a high density regime, µ q ≳ m π /2, has been observed by several lattice Monte Carlo simulations .
Until now, the two-point correlation function at nonzero µ q has been investigated [38,40,59,61].According to these works, it has been found that the two-point correlation function of a baryon has its imaginary-time reversibility broken once nonzero µ q is introduced.Thus, the dispersion relation is modified even in the hadronic phase, while we know the µ q -dependence analytically.This paper shows the µ-independence of some relevant parameters to scattering processes although the two-and four-point correlation functions clearly have µ-dependence.In the actual numerical calculation, it is hard to perform the exactly zero temperature simulation where we need infinite temporal lattice size.As we will see, however, the theoretical prediction at exactly zero temperature still holds in our numerical calculation at T ≈ 0.19T c .
Furthermore, using the breaking of the imaginary-time reversibility in the two-and four-point correlation functions at nonzero µ q , we propose a practical improvement method to obtain the more reliable value of the scattering phase shift.Thanks to the asymmetric behavior of the baryon correlation function, we can access long-τ data while avoiding the finite-T effect and can reduce a systematic error coming from inelastic contributions in the interaction potential.In the case of QC 2 D, the scattering parameters for pions and diquarks take the same values because of the Pauli-Gürsey symmetry, so that we can simultaneously obtain pion-pion scattering parameters, which are actually independent of µ in the hadronic phase as long as this improvement technique is employed.This paper is organized as follows.In Section 2, we formulate the (time-dependent) HAL QCD method in a way applicable to the hadronic phase with nonzero chemical potential and then clarify the µ-independence of hadron scatterings.In Section 3, we provide the HAL QCD method with a calculation strategy based on lattice QC 2 D Monte Carlo simulation.In Section 4, we therefrom demonstrate the µ-dependence as clarified in Section 2. In Section 5, we estimate a contribution from the inelastic channels to the interaction potential and phase shifts and then obtain our results for hadron scattering parameters in a long-τ regime by using asymmetric correlation functions for baryons at nonzero chemical potential.Section 6 is devoted to the summary and discussion.

Hadron scatterings at nonzero chemical potential
In this section, we present a theoretical discussion of hadron scatterings at T = 0 in QCD with nonzero µ and arbitrary N c .The system considered here stays in the hadronic phase at zero temperature, even though we introduce the number operator associated with nonzero µ into the QCD action1 .In Section 2.1, we derive the dispersion relation of a single hadron at nonzero µ.After briefly reviewing the HAL QCD method for describing hadron scatterings at µ = 0 in Section 2.2, we extend its formulation to the µ ̸ = 0 case in Section 2.3.Note that the conclusions obtained in this section hold true irrespective of the type of µ, namely, no matter whether it is the quark chemical potential (µ q ) or the baryon chemical potential (µ B ).

Two-point correlation functions and dispersion relation of a hadron at small chemical potential
To investigate scatterings, we first need to obtain the dispersion relation of a single particle.At zero chemical potential and zero temperature, a relativistic particle such as a hadron is known to have an energy E = p 2 + m 2 , where p and m are its momentum and mass, respectively.Here, we derive the dispersion relation of a hadron in the hadronic phase with nonzero chemical potential.Such a dispersion relation can be obtained from the time dependence of a (Euclidean) two-point correlation function in QCD defined by where Z = Tr e − 1 T ( Ĥ−µ N ) and T denote the partition function and temperature, respectively.Also, N denotes the operator of the particle number conjugate to µ: N is the quark number operator if µ = µ q , while N is the baryon number operator if µ = µ B .The operator Ô(τ ) is an arbitrary Heisenberg operator at Euclidean time τ .Here, we assume that Ô(τ ) In the path-integral formulation, the two-point correlation function is given by where O[U, ψ, ψ](τ ) is the corresponding function of field configurations to Ô(τ ) at τ .In the zero-temperature limit T → 0, it reads Using the Euclidean-time dependence of the Heisenberg operator at finite chemical potential and the commutation relation with N , we obtain (2.4) Thus, C(µ; τ ) can be described as where E n is the energy of the single hadronic state |n⟩ at µ = 0.When we set Ô(0) and Ô † (0) to the annihilation and creation operator of a single particle with momentum p, Eq. (2.5) at large τ gives the dispersion relation of a single hadron at small chemical potential as Thus, the meson energy is independent of µ since n O = 0, while the baryon energy depends on µ because of n O ̸ = 0. From Eq. (2.6), we observe an ambiguity arising from the definition of the "effective mass" of hadrons, in particular baryons, at nonzero µ.One definition is m eff := m; in this case, the effective mass is independent of µ, which in turn can be regarded as an energy shift, namely, E(p = 0, µ) = E(p = 0, 0) − µn O .This definition naturally matches the Silver Blaze phenomena.In fact, in the context of the lattice numerical studies on dense QC 2 D, Muroya et al. investigated the hadron mass spectrum using this definition [38].Another definition is m eff := E(p = 0, µ) or, equivalently, m eff = m − µn O .In this case, the exponential decay rate of the two-point correlation function in the imaginary time direction corresponds to the effective mass, which is in turn related to the pole mass.In this definition, it is possible that a baryon with n O > 0 becomes lighter than meson2 .
Several recent works [40,59,61] on dense QC 2 D have employed the effective mass as defined in the latter way.From now on, we call the former and latter definitions µ-independent scheme and µ-dependent scheme for the hadron effective mass, respectively.Note that the difference between the two schemes is just apparent; we can convert the former into the latter analytically and vice versa.
We must stress that the above discussion is applicable only when the particle (quark or baryon) number remains conserved; that is, the eigenvalue of N is zero ( N |0⟩ = 0).In the case of 3-colour QCD, it is widely believed that this is true for 0 < µ q < m N /3 ≈ µ c since µ q = m N /3 is a threshold value to create baryons dynamically.In the lattice QCD theory, on the other hand, it has been proven only for 0 < µ q < m π /2 at the T → 0 limit analytically [12][13][14].Actually, if we calculate the value of the quark number operator for the generated configuration using the lattice Monte Carlo method, then it takes almost zero for each configuration (|Ω i ⟩), i.e., ⟨Ω i |n q |Ω i ⟩ = 0 without taking the configuration average, in the range of 0 < µ q < m π /2, but in the m π /2 < µ q < m N /3 regime, each configuration gives a nonzero complex value of the quark number operator.⟨n q ⟩ = 0 is nevertheless suggested for m π /2 < µ q < m N /3 after taking the ensemble average [2,15], although it has yet to be proven.The lattice Monte Carlo calculations may well reproduce the physical quantity after taking such an ensemble average.Therefore, we conjecture that our formulation is valid also for m π /2 < µ q < m N /3 ≈ µ c and eventually in the whole µ q regime of the hadronic phase.In contrast to the hadronic phase, our formulation is invalid for µ q > µ c , since it is believed that a transition to the phase with ⟨n q ⟩ ̸ = 0, namely, the nuclear superfluid phase, occurs.
In the case of QC 2 D, for which we will demonstrate the simulation in the later part of this paper, the critical value µ c amounts to m π /2 since m N = m π at µ q = 0, so that the applicable range is 0 < µ q < m π /2.As for µ > µ c , the µ-dependence of the correlation function becomes more complex because of the spontaneous breaking of the particle number conservation, which changes the properties of the vacuum |0⟩ as referred to in Eq. (2.5).

Brief review of the HAL QCD method at zero chemical potential
Here, we give a brief review of the HAL QCD method, which has been used to investigate hadron scatterings at µ = 0 [18].The HAL QCD method is one of the techniques that allow us to extract S-matrix elements such as phase shifts via the interaction potential between two hadrons.In this work, we focus on the scattering only between the same hadrons and take the center-of-mass frame.
Let us start with the equal-time Nambu-Bethe-Salpeter (NBS) wave function defined as where O(x, τ ) denotes an operator located at (x, τ ), and |HH; W ⟩ is the two-body hadron state with energy W .Then, the energy satisfies the dispersion relation, The NBS wave function satisfies the Helmholtz equation as Another important property of the NBS wave function is that the asymptotic behavior of the radial part of the NBS wave function reads [17,[63][64][65]] for a given angular momentum l.Here, δ l (p) in the NBS wave function is equivalent to the phase shift of the S-matrix in the two-body scattering case.
The question is how to find the NBS wave function.We first identify the nonzero value of (p 2 + ∇ 2 )Ψ W (r) at finite r as an interaction potential term [16,17]: (2.10) The interaction potential U (r, r ′ ) can be derived from the so-called R-correlator, R(r, τ ) [18], which is composed of the four-point and two-point correlation functions.Here, the reduced mass is defined by m = m eff /2 with the effective mass (m eff ) of each hadron in the centerof-mass frame3 .Although the details will be explained in the next subsection, the basic strategy of the derivation is as follows.The R-correlator at large τ can be expressed by a linear combination of the NBS wave functions: where n labels the energy level, which is discretized at finite volume as shown just below.
Using this fact, we can obtain the interaction potential, U (r, r ′ ), from the following equation at large τ : (2.12) Then we can compute the NBS wave functions as solutions of the Schrödinger equation which contains the obtained potential U (r, r ′ ), and extract the scattering phase shift δ l (p) from their asymptotic behavior.
In general, to obtain the interaction potential U (r, r ′ ), we have to prepare an infinite number of the R-correlators with different combinations of the coefficients { Ãn } in Eq (2.11).In practice, we take the derivative expansion of U (r, r ′ ) and truncate the expansion to the leading order.Then the leading-order (LO) potential V LO (r) with U (r, r ′ ) ≃ V LO (r)δ (3) (r − r ′ ) can be extracted from a single R-correlator.In this case, the Schrödinger equation in the radial direction at a given orbital angular momentum l reads The obtained function, ψ W l (r), approximately satisfies the asymptotic behavior given by Eq. (2.9), so that we can obtain the approximated scattering phase shift.The phenomenologically relevant parameters such as the scattering length a H and the effective range r eff defined by the effective range expansion, can thus be computed with good accuracy, since the LO approximation to the interaction potential is valid in the low-energy regime.This is the basic idea of the HAL QCD method for obtaining the relevant parameters of the scattering processes in the infinite volume limit.
In finite volume, the momentum p is discretized to p n , so that the energy also takes the discrete value W n = 2 p 2 n + m 2 as introduced in Eq. (2.11).It is known that the phase shift δ l (p n ) at energy W n in finite volume is the same as the one in infinite volume [66][67][68].On the lattice, furthermore, we can calculate the R-correlator using the lattice simulations and then obtain the interaction potential V LO (r) nonperturbatively.The numerical data for V LO (r) are not continuous because of the finite lattice spacing, so that we estimate its continuous function form of r by interpolation.Finally, we solve the Schrödinger equation numerically and then compute the relevant parameters of the scattering processes.

HAL QCD method at small chemical potential
Now, we generalize the formulation of the time-dependent HAL QCD method to the case of finite µ.In this subsection, we confine ourselves to the finite spatial volume system in the T → 0 limit, although taking the infinite-volume limit does not change the discussion in this section.Thus, the momentum takes a discrete value p n .
Following Ref. [65], it is straightforward to conclude that the asymptotic behavior of the NBS wave function, Eq. (2.9), still holds in a regime of nonzero µ as long as the system stays in the hadronic phase (see Appendix A in detail).The µ-independence of scattering parameters essentially comes from this property.However, the HAL QCD potential, Eq. (2.10), can depend on µ via the reduced mass m, of which the µ-dependence we know analytically.It is still worth showing which quantities will be modified by the µ-insertion in the actual calculation processes shown in the previous subsection.It will also be helpful to see the validity of the present formulation.
We now consider the following quantity, namely, the R-correlator at finite µ, where C(µ; τ ) and F (µ; r, τ ) denote the two-point and four-point correlation functions, respectively.The four-point correlation function is defined by The operators O † (0)O † (0) represent the source of the two-body hadron state at τ = 0.The R-correlator is essentially the four-point correlation function, which is normalized by the two-point correlation functions to reduce statistical and systematic errors.Using the similar expansion for the two-point correlation function in Eq. (2.5), the four-point correlation function can be expanded as where n labels the energy level of the two-body hadron state |HH; W n ⟩, and Ψ Wn (r) = ⟨0| Ô(r, 0) Ô(0, 0) |HH; W n ⟩ denotes the NBS wave function with the state |HH; W n ⟩.Here, we put The abbreviation in the second line of the equation denotes the contribution from inelastic states, such as three-body states and two-body states composed of excited particles, which are suppressed at large τ .Via Eqs.(2.5) and (2.17), the R-correlator at large τ reads where with the single-hadron state with zero momentum |H; p = 0⟩, and Ãn = A n /C 2 .Note that the µ-dependence is the same between the numerator and denominator, so that the R-correlator is independent of µ4 .The right-hand side of Eq. (2.18) satisfies the Schrödinger equation as where Here, it is necessary to keep in mind that the reduced mass m = m eff /2 in the denominator of Eq. (2.19) can depend on µ if we adopt the µ-dependent scheme for the effective mass, while m in W n represents the effective mass at µ = 0.
Equation (2.19) can be rewritten as where is used.Therefore, R(µ; r, τ ) satisfies the following relation: for large τ where the inelastic scattering contributions are sufficiently suppressed.In the LO analysis, where we approximate the potential as the local one, U (r, r ′ ) ≃ V LO (r)δ (3) (r−r ′ ), the above equation reads Once we obtain the interaction potential, the later procedure for obtaining information on the scattering process is the same as the µ = 0 case: Solve the Schödinger equation (2.13) and then take the asymptotic limit of the resultant NBS wave function to obtain the phase shift and the other parameters.Now, we can see that the physical information on the scattering process does not depend on the choice of m, that is, the scheme for the effective mass at nonzero µ.The value of m alters the LO potential by a factor.According to Eq. (2.13), all terms have the same factor 1/ m, so that the solution of the differential equation, i.e., the NBS wave function, is independent of µ5 .This indicates that the phase shift of the S-matrix does not depend on µ for any scattering.
Note that following the same line of argument for deriving the µ-dependence of the dispersion relation, one can show the µ-independence of the scattering parameters unless a spontaneous breaking of the conservation of N occurs.Thus, after the phase transition to the broken phase, the R-correlator would inevitably begin to have nontrivial µ-dependence through the vacuum.
3 Calculation strategy in QC 2 D In this section, we perform the numerical simulation based on the theoretical method developed in the previous section.This is the first ab initio study of hadron scattering in QCD-like theory at nonzero µ q .Let us recall that the discussions in the previous section hold true for any number of colours of the SU(N c ) gauge theories.Here, we perform the lattice simulation in QC 2 D, which is one of the sign-problem-free QCD-like theories even at nonzero quark chemical potential.The one thing to keep in mind in the actual numerical simulation is that we take a finite temporal extent, corresponding to a nonzero temperature simulation.Our discussions in Section 2, which are presented at exactly zero temperature, are nevertheless worth examining in an actual calculation.
3.1 Flavour symmetry in 2-flavour QC 2 D Before going to lattice numerical calculations, we briefly give a summary of the flavour symmetry of our model.Here, we consider 2-flavour QC 2 D, where two quarks with different flavours have a degenerate mass.In the case of 3-colour QCD, the massless 2-flavour QCD has SU(2) L × SU(2) R × U(1) flavour symmetry.In the case of 2-colour QCD, however, it is enhanced to SU(4) symmetry because of the pseudo-reality of quarks.The extended flavour symmetry is called the Pauli-Gürsey symmetry [35,36].Once we introduce a small mass parameter, then Sp(4) ≃ SO(5) symmetry remains and five pseudo-Nambu-Goldston bosons appear, namely, three pions6 , a diquark, and an anti-diquark.Thus, these three mesons and two baryons have a degenerate mass at µ q = 0.
If we introduce a nonzero quark chemical potential, then the SO(5) symmetry breaks down to SU (2) V ×U (1) B .There still remains meson-baryon symmetry, but as we explained, the correlation functions for baryons (diquarks and antidiquarks) have µ q -dependence while the ones for mesons (pions) do not.Then, it is interesting to examine S-wave scatterings of two pions with isospin I = 2 (denoted as "ππ") and of two scalar diquarks (denoted as "DD") and compare them in each calculation step.
Beyond µ c , furthermore, the baryon symmetry, U (1) B , is spontaneously broken due to the diquark condensation, and then finally Sp(1) V ≃ SU (2) V remains in the regime of massive quarks and large µ q .In such a regime, our formulation cannot be applied since the vacuum |0⟩ is altered.

Lattice Setups
Let us explain our lattice setup for numerical simulations.We use the same lattice setup as that used in our previous works [48, 52, 55-57, 59, 60].The lattice gauge action used in this work is the Iwasaki gauge action [69], which is composed of the plaquette term with W 1×1 µν and the rectangular term with W 1×2 µν , with c 1 = −0.331and c 0 = 1 − 8c 1 .Here, β g = 4/g 2 0 in the 2-colour theory, and g 0 denotes the bare gauge coupling constant.
As for a lattice fermion action, we use the two-flavour Wilson fermion action, Here, the indices 1, 2 denote the flavour label, and µ q is the quark chemical potential.The Wilson-Dirac operator including the number operator, ∆(µ q ), is defined by where κ is the hopping parameter.Note that the γ 5 -Hermiticity is explicitly broken at nonzero µ q and that the Wilson-Dirac operator satisfies ∆ † (µ q ) = γ 5 ∆(−µ q )γ 5 .
We found that the superfluidity emerges at µ c /m 0 π ≈ 0.5 [52] as predicted by the chiral perturbation theory (ChPT) [70].It is natural to use µ q /m 0 π as a dimensionless parameter of density since the critical value µ c can be approximated as m 0 π /2 no matter how much the value of m 0 π is in a numerical simulation.Table 1 is the summary of the numbers of configurations and sources in the measurements of the two-point and four-point correlation functions.To measure the correlation functions, we take the wall source.Thus, the source operator at τ = τ 0 is defined by7 q w (τ 0 ) := y q(y, τ 0 ). (3.4) To improve the statistics, we take multiple timeslices of τ 0 at nonzero µ q .The statistical errors of our results are estimated by the jackknife method.

Two-point correlation functions for a pion and a diquark in QC 2 D
We calculate the two-point correlation functions of a pion and a scalar diquark given by where Here, K = Cγ 5 τ 2 is the anti-Hermitian operator which satisfies K 2 = −1, K † = −K, and C = γ 2 γ 4 is the charge conjugation operator in Euclidean spacetime.Note that in this notation, The Pauli matrix τ 2 acts on the colour index and hence is equivalent to iϵ ab .The index w in q w (τ 0 ) (q = u, d, ū, d) denotes the wall-source.
Taking all possible contractions in Eq. (3.5), the two-point function of a pion and of a diquark can be calculated by the following expression [40]: where Tr and subscript T denote the trace and the transpose for the colour and spinor indices, respectively.Here, S N denotes the (normal) quark propagator, In both the pion and diquark two-point functions, there is no disconnected diagram in the hadronic phase.The situation would be changed in the superfluid phase where scalar diquarks form a condensate.Then, we have to take the contributions from anomalous propagation for quark-to-quark or antiquark-to-antiquark into account.See Refs.[40,59] for the calculations of the two-point correlation functions in the superfluid phase.

Four-point correlation functions for pions and diquarks in QC 2 D
The four-point correlation functions for pions and diquarks are given by (3.9) Figures 2 and 3 present the all possible contractions for F π (r, τ ) and F D (r, τ ), respectively.As in the case of the two-point correlation functions, there is no disconnected diagram as long as we consider ππ and DD scatterings.Furthermore, there is no propagation between different space coordinates, (x + r) ↔ x, at the same time slice.Such a simple situation would not hold if we should perform a similar calculation in the superfluid phase.In the numerical calculations of the convolution (3.9), we utilize the fast Fourier transformation (FFT) library (see Appendix B for more details).Now, we can obtain the R-correlator numerically defined in Eq. (2.15).The next nontrivial calculation process is the approximation of the interaction potential U (r, r ′ ) to its LO potential V LO (r).To obtain the S-wave scattering amplitudes, we have to extract the LO potential from the S-wave component of the R-correlator 8 .The problem is that there is no longer continuous rotational symmetry on 3-dimensional spatial lattices.Instead of the continuous rotational symmetry, discrete rotational symmetry called the cubic group O holds on a lattice [71].The cubic group has an irreducible representation called A 1 representation, which is an alternative to the S-wave.We thus project the R-correlator onto the A 1 representation via where gr is the rotation of r by the element g, and χ A 1 (g) is the character of the A 1 representation.Finally, we can extract the LO potential V LO (r) from this projected Rcorrelator. 8The usual potential is independent of the angular momentum l.However, the LO potential in the HAL QCD method has an implicit l dependence.The reason is as follows.In the derivative expansion of the non-local potential, some higher-order terms include the angular-momentum operators L = −ir × ∇ because of the rotational symmetry.One example is V L 2 (r) L2 δ (3) (r − r ′ ).Once we choose specific l and its z component lz, the angular-momentum operators act trivially as (l(l + 1)V L 2 (r)δ (3) (r − r ′ ) for the above example).Then such higher-order terms become local and absorbed into the LO potential (see Section 4 in Ref. [17] for the case of the nucleon-nucleon potential).4 Simulation results

Two-point correlation function and effective mass
Figure 4 shows our numerical results for the two-point correlation function at different µ q .The pion two-point correlation function (left panel) is not only symmetric under time reversal, but also independent of µ q .On the other hand, the diquark two-point correlation function (right panel) becomes asymmetric at nonzero µ q and has its slope changed with increasing µ q .To obtain the effective mass, therefore, we use a single cosh function for the pion and a single exponential function for the diquark (antidiquark) as the fitting function.In the latter case, we define the antidiquark (diquark) propagator as the one in the forward (backward) temporal direction.
Recall that we have two schemes for the effective mass, namely, the µ-independent and µ-dependent schemes as discussed in Section 2.1.To extract the effective mass from each scheme, the fitting function in a long τ regime can be given as follows: µ-dependent scheme: Figure 5 depicts our results for the effective masses in the µ-dependent scheme.At µ q = 0, the masses of the diquark, the antidiquark, and the pion are degenerate, which comes from the unbroken part of the Pauli-Gürsey symmetry of 2-flavour QC 2 D [35,36].Our results for the pion and the diquark agree with m eff = m − µ q n q (black dashed lines): The pion mass is independent of µ q while the diquark mass behaves as m − 2µ q .As for the antidiquark, the results slightly deviate from the linear behavior at large µ q , a feature caused by the effect of the periodicity on finite lattices.Thus, as evident from the right panel of Figure 4, it is getting more difficult to secure the data for the antidiquark correlation function at sufficiently long τ as µ q increases.From now on, we focus only on the pion and diquark, which are the lightest meson and baryon in QC 2 D. The agreement of our numerical data with the lines of m eff = m − µ q n q indicates that the pion and diquark effective masses in the µ-independent scheme are actually independent of µ q .
Figure 5. Effective masses at each µ q for a pion (red circles), a diquark (magenta squares), and an antidiquark (brown crosses).The black dashed lines correspond to m 0 π ± 2µ q .

R-correlator at nonzero chemical potential
Our next task is to examine the R-correlator, which is the essential quantity in the timedependent HAL QCD method.We define the ππ and DD R-correlators as where F π (µ q ; r, τ ) (F D (µ q ; r, τ )) is the four-point correlation functions of pions (diquarks) defined in Eq. (3.9).We should emphasize that, as in the case of the diquark two-point function, the DD two-body states contribute to the diquark four-point correlation function in the backward temporal direction.Thus the DD R-correlator is composed of the diquark two-and four-point correlation functions in the backward direction, C D (µ q ; aN τ − τ ) and F D (µ q ; r, aN τ − τ ), respectively.The blue circles, green crosses, and red squares correspond to the results at µ q /m 0 π = 0.00, 0.16 and 0.32, respectively.To see the plots easily, the data are truncated at the timeslice where the R-correlator is maximally affected by the finite-T effect: τ /a = 16 for all µ q for the ππ system, and τ /a = 16, 21, and 26 for µ q /m 0 π = 0.00, 0.16 and 0.32, respectively, for the DD system.
Figure 6 depicts the logarithmic plot of the ππ and DD R-correlator summed over the spatial coordinates r R(µ q ; r, τ ) at each timeslice.We find that at every timeslice except τ /a = 15 and 16, the results for the ππ system with different µ q agree with each other, which is obvious because of the µ q -independence for the pion four-and two-point correlation functions predicted in our formulation.At τ /a = 15 and 16, the deviation of the data for µ q /m 0 π = 0.32 from those for µ q /m 0 π = 0.00, 0.16 is attributable to statistical fluctuations.Furthermore, we can see the rising behavior of the ππ R-correlator in τ /a ≳ 11, in which the data suffer from the finite-T effect and thus cannot be used for our analyses.
For the DD system, we can also see the agreement between the results for different µ q until τ /a ≈ 10.In contrast to the ππ case, however, the rising behavior appears later for larger µ q in this system: τ /a ∼ 11, 17, and 22 for µ q /m 0 π = 0, 0.16, and 0.32, respectively.This indicates that the finite-T effect on the DD R-correlator is suppressed by nonzero µ q .Such suppression can be explained as follows.The finite-T effect on the DD R-correlator comes from the contribution of the diquark two-point and four-point correlation functions in the forward temporal direction, which can be expanded as Eqs.(2.5) and (2.17) with the one-body and two-body antidiquark states, respectively.The exponent of these terms can be expressed as W − n Dµ q , where W is the energy at µ q = 0 and n D = −2 (n D = −4) for the two-point (four-point) correlation function.For a large µ q , therefore, such terms are suppressed at earlier τ in the forward direction; in other words, these contributions appear at later τ in the backward direction.The details of the τ -dependence will be revisited in Section 5.
From now on, we focus on the µ q -dependence of the potentials and phase shifts that can be obtained from the R-correlators at τ /a = 8.Here, we take this timeslice to suppress the contributions coming both from the excited states and the periodicity at µ q = 0 as much as possible.Figure 7 exhibits the ππ (left panel) and DD (right panel) R-correlators for different µ q .We find that both results for the ππ and DD systems have no µ q -dependence.Thus, the µ q -dependence of the DD four-point correlation function cancels with the one of the diquark two-point correlation function squared as expected from Eq. (2.18).

Leading-order potential and scattering phase shift
Finally, let us investigate the shape of the interacting potential, Eq. (2.23), and the phenomenological parameters for scattering processes at different µ q .Here, we have to be careful about the scheme dependence of the effective mass, which has relevance to the choice of m in Eq. (2.23).If we take the µ-independent scheme, i.e., m = m/2 in Eq. (2.23), then it is obvious that the LO potential does not depend on µ q even for DD scattering since, as we have already shown, our R-correlator is independent of µ q .Therefore, here, we take the µ q -dependent scheme, i.e., m = E(0, µ q ), and see the potential shape at different µ q .Furthermore, τ -dependence of the potential shape will be discussed in Section 5.
In Figure 8, we present the leading-order potentials, V LO (r).First, we find that the ππ and DD potentials are the same at µ q = 0, which can be explained by the unbroken part of the Pauli-Gürsey symmetry.For the ππ system shown in the left panel, the three potentials with different values of µ q agree within the error bars as they should.On the other hand, the DD potential clearly depends on µ q : As shown in the right panel of Figure 8, the repulsive  core around the origin and the attractive pocket at r ≈ 1.2 fm are getting stronger and deeper, respectively, with increasing µ q .Furthermore, in the long r regime, the results show an appreciable discrepancy between the potential at µ q /m 0 π = 0.32 and the others, which is not necessarily predictable.We have nevertheless found that the difference comes from the decrease of the reduced mass m at large µ q , which acts to scale up the potential (2.23).Indeed, the potentials multiplied by the reduced mass m(µ q )V LO (r) agree with each other as shown in Figure 9. Thus, we can conclude that the scheme dependence of the effective mass on the interaction potential appears only in the overall factor.
To see the phase shift of the S-matrix, we first perform the fitting for the obtained potential data using the following fitting function, V (r) = a 0 e −(r/a 1 ) 2 + a 2 e −(r/a 3 ) 2 + a 4 e −(r/a 5 ) 2 + a 6 e −(r/a 7 ) 2 .
(4.5)Then, we solve the Schrödinger equations in the radial direction, Eq. (2.13), where the angular momentum is set to l = 0 for the S-wave channel.The resultant radial part of the NBS wave function leads to the phase shift δ(p).In Figure 10, we show our results for p cot δ(p).For both the ππ and DD systems, the results for different µ q agree with each other within the error bars, which is consistent with the prediction in Section 2.3.It is also found that the ππ and DD phase shifts match each other within the error bars.This is reasonable because, at µ q = 0, the pion and the scalar diquark belong to the same multiplet of the unbroken part of the Pauli-Gürsey symmetry, which remains even at µ q ̸ = 0 as long as the system stays in the hadronic phase.Thus, all two-body scatterings for pions, diquarks, and anti-diquarks have the same value of the phase shift in the hadronic phase at nonzero µ q .
5 Analysis of leading-order DD (ππ) potential in a long-Euclidean-time regime at nonzero chemical potential As seen in the right panel of Figure 6, we have found that the finite-T effect on the DD R-correlator gets suppressed by nonzero µ q .In this section, we obtain a more reliable value of the phase shift by making the most of the property at nonzero µ q .Thus, we examine the LO potential (2.23) for the DD system that can be obtained from the R-correlator in a long-τ regime at the largest value of µ q and derive the p cot δ(p) from the resulting NBS wave function.We have shown that the obtained value of p cot δ(p) is the same for both ππ and DD (also D D) in the hadronic phase.By taking the long-τ data at nonzero µ q , we can reduce systematic errors coming from the inelastic contribution in our formula (2.17) and also finite-T effects.Throughout this section, we take the µ-independent scheme of the reduced mass in our calculations.Figure 11 shows the LO potential at τ /a = 8, 11, 14, and 17, in which region the finite-T effect is under control.Note that the results at τ /a = 8 are the same as those in both panels of Figure 8 for µ q = 0.For τ /a = 8, the obtained potential shape suggests a repulsive core at short distances with a rather small attractive pocket.On the other hand, the shape of the potential for τ /a ≥ 11 changes from that for τ /a = 8: The repulsive core gets slightly larger, while the attractive pocket disappears.Furthermore, no drastic difference can be seen among the three potentials for τ /a ≥ 11.We expect that the contribution from the inelastic states causes such discrepancy between the potentials at different τ .Although we assume the inelastic contribution to be negligible in deriving the LO potential in Section 4.3, the τ dependence of the LO potential suggests that the resultant systematic errors still remain even at τ /a = 8.
The repulsion at short distances without any attractive pockets in the LO potential at later τ (e.g., τ /a = 17) is qualitatively consistent with the earlier results of quenched QC 2 D simulation, namely, QC 2 D simulation without dynamical quarks [72].Furthermore, such a repulsive core has also been seen from the results for the I = 2 pion-pion potential obtained for 3-colour QCD with µ q = 0 [73,74].A good observable to see a precise τ -dependence is p cot δ(p).We extract this quantity by fitting the LO potentials and solving the Schrödinger equations in the radial direction, as explained in Section 4.3.The details of the fitting results are presented in Appendix.C. The results for p cot δ(p) are shown in Figure 12.Now, we can observe that the results at τ /a = 14 and τ /a = 17 are consistent with each other, but the result at τ /a ≤ 11 deviates from them beyond the statistical error bars.Furthermore, the scattering length a H and effective range r eff obtained by fitting the resultant p cot δ(p) to Eq. (2.14) are listed in Table 2.These results indicate that the contribution from the inelastic states would be sufficiently suppressed if we can take τ /a ≳ 14.

Summary and discussion
In this paper, we examined the chemical potential dependence of hadron scatterings.We focus on the system with small µ, which still stays in the hadronic phase where the baryon number has yet to be violated.As is analytically known, the shape of the hadron twopoint correlation functions depends on µ even in the hadronic phase, and accordingly the dispersion relation of a single hadron is modified.Given such a µ dependence, it might be possible to see that the hadron mass spectra are unchanged, a property that is consistent with the Silver Blaze phenomena.In this paper, we also develop a theoretical discussion of what are µ-(in)dependent quantities by extending the time-dependent HAL QCD method at µ = 0 to µ ̸ = 0.
We first show the dispersion relation of a single hadron at small µ from the Euclideantime (τ ) dependence of its two-point correlation function, which gives E(p, µ) = p 2 + m 2 − µn O with the corresponding quantum number of the hadron (operator) n O .We can consider two schemes for the effective mass: the µ-independent scheme m eff = m and the µ-independent scheme m eff = m − µn O .
Next, we formulated the time-dependent HAL QCD method for describing hadronhadron scatterings at small µ, where we utilize the above effective mass in the reduced mass ( m).We found that the R-correlator, composed of the four-point and two-point correlation functions, does not depend on µ because the µ-dependence of the two-point and four-point correlation functions is canceled out.Then, the HAL QCD potential, which can be calculated from the R-correlator, depends on µ only through the reduced mass.Finally, we show that the asymptotic behavior of the NBS wave function and the scattering amplitude are independent of µ in the hadronic phase.Thus, the physical information on scattering processes does not depend on the choice of the reduced mass, that is, the scheme for the effective mass.
To demonstrate the above analytical predictions, we analyzed the S-wave scatterings of two pions with isospin I = 2 (ππ) and two scalar diquarks (DD) by applying the HAL QCD method to QC 2 D at nonzero quark chemical potential µ q in the µ-dependent scheme.This is the first study on the numerical calculations for the hadron scattering in QCD-like theory at finite µ q .The R-correlators for different µ q match each other within the error bars for both the ππ and DD systems unless the finite-T effect is comparable.We also have found that the DD R-correlator for larger µ q suffers less from the finite-T effect thanks to the asymmetric properties of the correlation functions.Furthermore, our results indicate that the LO ππ potential is independent of µ q while the LO DD potential depends on µ q through the factor 1/ m(µ q ).Finally, the results for the scattering phase shifts do not depend on µ q .Furthermore, the values of the relevant parameters of scattering for both systems agree with each other because of the partial Pauli-Gürsey symmetry of 2-flavour QC 2 D. The behavior is consistent with the predictions from our formulation.
In addition, we have succeeded in obtaining the reliable DD potential and phase shift by using the DD R-correlator in a long-τ regime at µ q /m 0 π ≈ 0.32.Due to the mesonbaryon symmetry of 2-flavour QC 2 D and the µ-independence of the phase shift, the results are the same for both the DD and the ππ systems in the whole µ q regime in the hadronic phase.The resultant potential shape is characterized by the repulsive potential without any attractive pocket, which is qualitatively consistent with the previous results of the quenched QC 2 D and 3-colour QCD with the heavy pion mass.The results for the scattering phase shifts indicate that the contribution from the inelastic states is sufficiently suppressed in τ /a ≳ 14.The precise determination of the scattering parameters in this system using a larger temporal volume simulation combined with a higher-order analysis in the timedependent HAL QCD method is a future work.
The suppression of the finite-T effect by the nonzero-µ insertion, as observed in the right panel of Figure 6, occurs for the hadron operators with a positive quantum number for any QCD-like theories.Once we obtain the correlation functions at µ ̸ = 0 and know their µ-dependence analytically by following the derivation in Section 2.1 or 2.3, we can estimate the ones at µ = 0 in a manner to keep the finite-T effect suppressed.For example, if we consider 3-colour QCD with nonzero isospin chemical potential µ I [75], we can likewise reduce the finite-T effect, that is, the periodicity effect of the correlation functions of π − mesons.The application to other systems and/or other types of chemical potentials is left for our interesting future works.
The discussion in this paper is valid unless the spontaneous breaking of the conservation occurs.We expect that after the phase transition from the hadronic to superfluid phase, the R-correlator would have nontrivial µ dependence through the vacuum.The examination of such behavior is one of the important plans for our future studies.

C Details of fitting potentials
In this appendix, we show the details of fitting potentials presented in Section 5, such as the setup and the results.As a fit function, we utilize the sum of four Gaussians given by V (r) = a 0 e −(r/a 1 ) 2 + a 2 e −(r/a 3 ) 2 + a 4 e −(r/a 5 ) 2 + a 6 e −(r/a 7 ) 2 , (C.1) where we assume that a 1 < a 3 < a 5 < a 7 .We employ the uncorrelated fit in this analysis.17 (lower right) for µ q /m 0 π = 0.32 in the µ-independent scheme, where we can regard these results as those for ππ potentials at µ q = 0.
The results for the best-fit parameters in Section 5 are listed in Table 3.Also, Figure 13 depicts the comparison of the fitted results and the raw data for V LO (r).We observe that all fit results except for τ /a = 14 and 17 agree with the original data, which implies that the fitting works well.
On the other hand, the original data and the fit results for τ /a = 14 and 17 deviate from each other in the long-range region: The original data in r ≳ 3 fm do not seem to have reached zero while the fit results converge to zero.We consider that such deviation comes from the statistical fluctuations around zero.Indeed, we find that DD potential fluctuates around zero at long distances when the Euclidean time is varied around those timeslices.Therefore, we here neglect such a long-range deviation.

Figure 1 .
Figure 1.Conjectured QCD phase diagram.The horizontal axis denotes the quark chemical potential (µ q ) and n O = −N c for the meson, baryon, and antibaryon operators, respectively.On the other hand, if µ = µ B , then n O = 0, n O = 1, and n O = −1 for each of them.

Figure 6 .
Figure 6.ππ (left) and DD (right) R-correlator summed over the spatial coordinates r at each timeslice τ .The blue circles, green crosses, and red squares correspond to the results at µ q /m 0 π = 0.00, 0.16 and 0.32, respectively.To see the plots easily, the data are truncated at the timeslice where the R-correlator is maximally affected by the finite-T effect: τ /a = 16 for all µ q for the ππ system, and τ /a = 16, 21, and 26 for µ q /m 0 π = 0.00, 0.16 and 0.32, respectively, for the DD system.

Figure 7 .
Figure 7. ππ (left) and DD (right) R-correlators at τ /a = 8.The blue circles, green crosses, and red squares correspond to the results at µ q /m 0 π = 0.00, 0.16 and 0.32, respectively.Note that the green and red plots are slightly shifted horizontally for visualization, although they are still almost invisible in the figure because all the plots including the blue one are the numerical results for the same quantity.

Figure 8 .
Figure 8. Leading-order ππ (left) and DD (right) potentials at τ /a = 8.Note that the green and red points are almost invisible in the left figure because all the points including the blue ones are the numerical results for the same quantity.

Figure 9 .
Figure 9. Leading-order DD potentials at τ /a = 8 multiplied by the reduced mass.Note that green and red points are almost invisible because all the points including the blue ones are the numerical results for the same quantity.

Table 2 .
DD scattering length a H and effective range r eff derived using the LO potentials at τ /a = 8, 11, 14, and 17 for µ q /m 0 π ≈ 0.32.