Detecting cosmological phase transitions with Taiji: sensitivity analysis and parameter estimation

The direct detection of gravitational waves (GWs) by the LIGO and Virgo collaborations [1] has initiated a new era in observational astronomy, providing unprecedented access to astrophysical phenomena that remain invisible to electromagnetic observations. While ground-based detectors operate approximately between 10 Hz and 1 kHz, space-based interferometers will explore the milli-Hertz frequency band, where signals from various cosmological sources are expected to be present [2, 3].

One of the potential targets for space-based GW observatories is stochastic GW background (SGWB) produced by first-order phase transitions (FOPTs) in the early universe [4, 5]. These transitions occur when a system transitions discontinuously between different vacuum states separated by an energy barrier, resulting in the nucleation and expansion of bubbles of the new phase within the old phase [6−20]. In the Standard Model of particle physics, the electroweak phase transition is a crossover-type; however, many well-motivated extensions predict a first-order electroweak phase transition occurring at temperatures of $ T_* \sim 100 $ GeV [21, 22]. Such phase transitions could explain the observed baryon asymmetry of the universe through electroweak baryogenesis [23, 24], and might be connected to the production mechanisms of dark matter [25].

The Chinese space-based GW observatory Taiji [26, 27] is one of several proposed missions designed to detect GWs in the milli-Hertz frequency range. Similar to the European Space Agency's Laser Interferometer Space Antenna (LISA) [28], Taiji will consist of three spacecraft in a triangular formation, but with arm lengths of $ 3 \times 10^6 $ km compared to LISA's $ 2.5 \times 10^6 $ km. The Taiji constellation will follow a heliocentric orbit about $ 20^\circ $ ahead of Earth. Another Chinese space-based detector, TianQin [29], is designed with shorter arm lengths of $ \sim 10^5 $ km and will be in the Earth's orbit, providing complementary sensitivity in partially overlapping frequency bands.

Detecting the SGWB from FOPTs requires distinguishing this cosmological signal from foreground sources, primarily from galactic and extragalactic compact binary (ECB) systems. The unresolved population of double white dwarf (DWD) binaries in our galaxy forms a significant confusion foreground [30, 31], whereas the superposition of signals from ECB coalescences contributes an additional stochastic background [32, 33]. The Taiji mission, with its specific noise characteristics and orbital configuration, presents unique capabilities and challenges for separating these components.

The detection of an SGWB from FOPTs faces significant challenges, primarily due to the contamination of SGWBs from unresolved galactic compact binaries, particularly DWD systems [30, 34], and from extragalactic compact binaries [35]. The astrophysical SGWBs are sufficiently strong and become foregrounds acting as additional ``confusion noise" when conducting the detections of other GW signals in the same frequency band [36]. Therefore, the foregrounds must be carefully modeled and subtracted to reveal primordial signatures [37, 38].

Previous studies have investigated LISA's capabilities to detect SGWBs from FOPTs [39−46]. More recently, attention has turned to the complementary capabilities of Taiji and the potential for joint observations with LISA or TinQin [47−54]. However, a comprehensive analysis of Taiji's sensitivity to FOPTs, considering the latest noise models and foreground estimates, remains to be conducted.

In this study, we comprehensively assess Taiji's capability to detect SGWBs from FOPTs. Our analysis incorporates detailed modeling of the Taiji noise spectrum, including both instrumental noise and astrophysical foreground contributions from galactic DWD binaries and ECB systems. We implement a Bayesian framework to systematically explore the detectability of phase transition signals across a range of amplitudes and peak frequencies, determining the regions of parameter space where Taiji can make robust detections and provide precise parameter estimates. The remainder of the paper is organized as follows: In Section II, we present our models for the GW signal from FOPTs, the Taiji detector sensitivity, and the relevant astrophysical foregrounds. In Section III, we describe our Bayesian methodology and simulation framework. Finally, in Section IV, we summarize our findings and discuss their implications for probing beyond-Standard-Model physics with future space-based GW observatories.

In this section, we describe the key components of our analysis framework. We first present our model for the GW signal from FOPTs, followed by a detailed characterization of the Taiji detector's noise properties. We then discuss the two primary astrophysical foregrounds that impact the detection of cosmological signals: the galactic DWD confusion noise and the ECB background.

For our analysis of FOPTs as sources of an SGWB, we employ a simplified broken power-law spectral model derived from fitting to numerical simulations [11]. The GW energy density is

where the spectral shape function takes the form of

Here, $ \Omega_{\text{PT}} $ represents the peak amplitude of the SGWB, and $ f_{\text{PT}} $ is the peak frequency [11]. The GW power spectrum relates to the power spectral density at the detector through

where $ H_0 = 67.4 \, \text{km} \, \text{s}^{-1} \, \text{Mpc}^{-1} $ is the Hubble parameter today [55]. The peak frequency depends on the physical parameters of the phase transition:

where $ T_* $ is the temperature at which the phase transition occurs, $ H_* $ is the Hubble rate at that time, and $ R_* $ is the mean bubble separation.

The $ f^3 $ low-frequency behavior of $ {\cal{P}}(f) $ in Eq. (2) is characteristic of phase transitions with mean bubble spacing on the order of the Hubble radius, which produce the strongest signals [19, 20]. The high-frequency $ f^{-4} $ behavior approximates the falloff seen in numerical simulations near the peak [11]. This model captures the essential features of FOPT signals while reducing the parameter space to two physically meaningful parameters: $ \Omega_{\text{PT}} $ and $ f_{\text{PT}} $. For phase transitions in the temperature range of 100 GeV to 1 TeV (including the electroweak scale and many BSM scenarios), we expect peak frequencies between $ 10^{-4} $Hz and $ 10^{-2} $ Hz with peak amplitudes in the range $ 10^{-14} < \Omega_{\text{PT}} < 10^{-9} $ [41]. These signals fall squarely within Taiji's sensitivity band, making Taiji a promising detector for probing BSM physics through GWs from FOPTs.

The Taiji space-based GW observatory features three spacecraft in a triangular configuration with 3 million kilometer arm lengths, longer than LISA's 2.5 million kilometers [26, 48]. To extract GW signals from the raw measurements, Taiji employs sophisticated signal processing techniques known as time delay interferometry (TDI) [56, 57]. For our analysis, we focus on the interferometric data streams designated as the X, Y, and Z TDI variables, which represent combinations of phase measurements that substantially reduce laser frequency noise. We adopt several simplifications in our noise modeling approach: 1) we assume that the SGWB signal and instrumental noise are uncorrelated, 2) we model the noise as consisting of two primary components, and 3) we treat all spacecraft as identical with equal arm lengths forming an equilateral triangle with $ L = 3 \times 10^{6} $ km [58].

The two dominant noise contributions in the Taiji detector can be characterized by their power spectral densities (PSDs). The first component arises from the optical measurement system (OMS), which dominates at higher frequencies (see e.g., [59])

where $ P = 8 $ [60]. The second noise component comes from acceleration noise affecting the test masses, which dominates at lower frequencies.

where $ A = 3 $ characterizes the acceleration noise level [60].

With these noise components defined, we can express the noise auto-correlation in the X, Y, and Z channels as

where c is the speed of light, and $ f_* = c / (2 \pi L) $ defines a characteristic frequency of the detector geometry. The cross-correlation between different channels (e.g., between X and Y) is given by

For analytical convenience, we transform the X, Y, and Z channels into an alternative basis consisting of the channels A, E, and T.

This transformation is advantageous because it produces noise-orthogonal channels $ \mathrm{A} $ and $ \mathrm{E} $ with identical noise properties, whereas $ \mathrm{T} $ functions as a “null channel” with reduced sensitivity to GWs [61]. The noise power spectra in these channels can be derived as

and

To facilitate comparison with astrophysical and cosmological GW signals, we convert the noise spectral densities to equivalent energy spectral densities as

where $ \alpha \in \{\mathrm{A, E, T}\} $ denotes the channel, and $ H_0 $ is the Hubble constant. The noise spectral densities $ S_\alpha(f) $ for each channel are defined as

Here, $ {\cal{R}}_\alpha $ corresponds to the response function for the respective channel α. For this analysis, we employ the analytical expressions for these response functions as derived in [62].

The Milky Way hosts a vast population of DWD binaries, with population synthesis models suggesting that there are approximately $ 10^7-10^8 $ such systems throughout our galaxy [63, 64]. These binaries generate gravitational radiation primarily in the frequency band spanning from $ 10^{-5} $ to $ 10^{-1} $ Hz [65], which overlaps significantly with Taiji's detection window.

While Taiji will resolve individual signals from the strongest and closest sources, the vast majority of these binaries produce signals below the detection threshold. These unresolved systems generate a collective SGWB that manifests as an additional noise component in the detector, commonly referred to as the ''confusion noise'' or ''galactic foreground'' [66]. For a 4-year observation period, we approximate this galactic background using a broken power-law model

The parameters that best fit the detailed population model are $ A_1 = 3.98 \times 10^{-16} $, $ A_2 = 4.79 \times 10^{-7} $, $ \alpha_1 = -5.7 $, and $ \alpha_2 = -6.2 $ [67, 68]. This functional form captures the essential spectral features of the DWD background, particularly the high-frequency steepening that occurs as the number of contributing binaries decreases. This spectral break arises from physical constraints on binary orbital separations, which cannot be smaller than the combined radii of the component white dwarfs. The corresponding energy density spectrum, normalized to the critical density of the universe, is given by

Beyond our galaxy, the universe contains innumerable compact binary systems that collectively generate an SGWB. This cosmological signal differs fundamentally from the galactic foreground, as it represents the superposition of unresolved binary black hole and neutron star systems that have been distributed throughout cosmic history [69].

While current ground-based interferometers have not yet reached the sensitivity required to detect this background, space-based detectors operating at lower frequencies will probe a different portion of its spectrum. The ECB background is important for understanding the integrated merger history across cosmic time.

For our sensitivity analysis, we model this background with a characteristic power-law frequency dependence

This spectral shape emerges naturally from the inspiral phase of compact binaries, with the $ 2/3 $ power-law index reflecting the frequency evolution of binary systems dominated by gravitational radiation. We adopt an amplitude of $ A_{\mathrm{ECB}} = 1.8 \times 10^{-9} $ at the reference frequency $ f_{\mathrm{ref}} = 25 $ Hz [69].

Unlike the galactic foreground, this background exhibits no spectral breaks within the Taiji frequency band, as the contributing sources span a much broader range of masses, redshifts, and formation channels.

This section outlines our computational approach for evaluating the Taiji mission's capability to detect SGWBs from cosmological FOPTs following [70, 71]. Our numerical framework simulates Taiji observations spanning the full 4-year mission duration, incorporating with realistic duty cycle considerations (assuming 75% efficiency [70, 72, 73]), yields an effective 3-year observations. We segment the TDI measurements into $ N_c = 94 $ chunks of 11.5 days each [70, 71]. The frequency domain extends from $ 3 \times 10^{-5} $Hz to $ 0.5 $ Hz with approximately $ 5 \times 10^7 $ total data points at $ 10^{-6} $ Hz resolution.

For computational implementation, we transform the time-domain signal into frequency space:

Under the assumption of stationarity for both signal and noise components, the Fourier coefficients exhibit the following statistical properties:

The simulation generates synthetic observations by drawing complex Fourier coefficients from Gaussian distributions characterized by the appropriate power spectral densities. Specifically, at each frequency point, we construct:

Here, $ G_{ij}(M,\sigma) $ represents random samples from a Gaussian distribution with mean M and standard deviation σ. The total power at each frequency combines signal and noise contributions: $ D_i = S_i + N_i $. To account for statistical fluctuations, we generate $ N_\mathrm{c} $ independent realizations $ \{D_{i1}, D_{i2}, \ldots, D_{iN_\mathrm{c}}\} $ at each frequency and compute their ensemble average $ \bar{D}_i $. Figure 1 illustrates a representative simulated dataset, with injection parameters documented in Table 1.

Figure 1.  (color online) Frequency-domain representation of synthetic Taiji A-channel observations (blue). We also show the galactic DWD confusion noise (orange), contribution from ECB (green), and cosmological background from the FOPT (red) with parameters $\Omega_{\mathrm{PT}}=3.9\times 10^{-11}$ and peak frequency $f_{\mathrm{PT}}= 7\times 10^{-3} \mathrm{Hz}$. For reference, the Taiji detector's sensitivity is plotted in terms of $\Omega_{\mathrm{GW}}(f)$ as a gray curve.

Parameter	Prior	Injected value	Recovered value
A	$ {\cal{U}}(2.95, 3.05) $	$ 3 $	$ 3.002^{+0.007}_{-0.007} $
P	$ {\cal{U}}(7.99, 8.01) $	$ 8 $	$ 7.9990^{+0.0013}_{-0.0014} $
$ \log_{10} A_1 $	$ {\cal{U}}(-16, -15) $	$ -15.4 $	$ -15.39^{+0.04}_{-0.05} $
$ \alpha_1 $	$ {\cal{U}}(-6, -5.5) $	$ -5.7 $	$ -5.69^{+0.05}_{-0.05} $
$ \log_{10} A_2 $	$ {\cal{U}}(-6.5, -6) $	$ -6.32 $	$ -6.31^{+0.04}_{-0.04} $
$ \alpha_2 $	$ {\cal{U}}(-6.5, -6) $	$ -6.2 $	$ -6.19^{+0.04}_{-0.04} $
$ \log_{10} A_{\mathrm{ECB}} $	$ {\cal{U}}(-9, -8.5) $	$ -8.74 $	$ -8.69^{+0.15}_{-0.18} $
$ \alpha_{\mathrm{ECB}} $	$ {\cal{U}}(0.5, 1) $	$ 2/3 $	$ 0.68^{+0.04}_{-0.05} $
$ \log_{10} \Omega_{\mathrm{PT}} $	$ {\cal{U}}(-10.609, -10.209) $	$ -10.409 $	$ -10.408^{+0.004}_{-0.004} $
$ \log_{10} (f_{\mathrm{PT}}/\mathrm{Hz}) $	$ {\cal{U}}(-2.355, -1.955) $	$ -2.155 $	$ -2.154^{+0.002}_{-0.001} $

Table 1.  Summary of Bayesian analysis results for all model parameters. The table displays the uniform prior ranges ($ {\cal{U}} $) employed in our MCMC sampling, alongside the true parameter values used in synthetic data generation. The rightmost column presents the posterior estimates, showing median values with corresponding 90% credible intervals. The prior ranges are chosen to balance computational efficiency with statistical robustness while focusing on the theoretically motivated parameter space for FOPTs detectable by Taiji.

To enhance computational efficiency while preserving information content, we implement adaptive frequency binning. For frequencies below $ 10^{-3} $ Hz, we maintain the original resolution, whereas frequencies between $ 10^{-3} $ Hz and $ 0.5 $Hz are rebinned into 1000 logarithmically spaced intervals. This optimization reduces the dataset to 1971 frequency bins per segment. The rebinned data are calculated as

where the optimal weights are

Here, $ {\cal{D}}^\mathrm{th}(f_j, \vec{\theta}, \vec{n}) \equiv \Omega_{\mathrm{GW}} (\vec{\theta}, f_j) + \Omega_{\alpha} (\vec{n}, f_i) $ represents the theoretical model for the total energy density, which is an estimate of the variance of the segment-averaged data $ \bar{D}_j $ [71]. The parameter $ \vec{n} \equiv \{{\cal{A}}, P\} $ denotes the instrumental noise parameters, whereas $ \vec{\theta} \equiv \{A_1, \alpha_1, A_2, \alpha_2, A_\mathrm{ECB}, \alpha_\mathrm{ECB}, \Omega_\mathrm{PT}, f_\mathrm{PT}\} $ encompasses all astrophysical and cosmological signal parameters, including the galactic DWD foreground, ECB background, and the FOPT signal of interest.

We now provide a brief derivation of Eq. (24). Each data point $ \bar{D}_j $ has variance $ \text{Var}(\bar{D}_j) = D_{\text{th}}(f_j, \vec{\theta}, \vec{n}) $. For the binned estimator $ D_{(k)} = \sum\nolimits_{j \in \text{bin}\,k} w_j \bar{D}_j $, assuming uncorrelated data points within each bin, the variance is

To minimize this variance subject to the normalization constraint

we use the method of Lagrange multipliers. The Lagrangian is

Taking the derivative with respect to $ w_i $ and setting it to zero yields

Applying the normalization constraint in Eq. (26), we obtain

Substituting Eq. (29) back yields the optimal weights in Eq. (24).

Our statistical analysis employs a hybrid likelihood function combining Gaussian and log-normal components [71], namely,

The Gaussian part is

while the log-normal part is

The inclusion of the log-normal component in our likelihood function is crucial for properly handling the statistical properties of power spectral densities. When analyzing SGWB signals, the power spectral densities follow an $ \chi^2 $ distribution rather than a Gaussian distribution. Using a Gaussian likelihood exclusively in such cases can lead to biased parameter estimation, particularly for weak signals with a low signal-to-noise ratio [71]. The log-normal term better captures the right-skewed nature of the $ \chi^2 $ distribution while maintaining computational tractability. This hybrid likelihood approach has been widely adopted and validated in the literature for SGWB analyses (see e.g., [38, 71, 74]).

To quantitatively assess the detectability of phase transition signals, we employ two complementary model selection metrics: the Bayes factor (BF) and the Deviance Information Criterion (DIC). The Bayes factor represents the ratio of evidence between competing models, providing a direct measure of the relative probability of the models. Specifically, we define BF as

where $ {\cal{Z}}_\text{FOPT} $ is the evidence for the model including a PT component, and $ {\cal{Z}}_\text{null} $ represents the model with only astrophysical foregrounds and instrumental noise. The values of $ \ln (\mathrm{BF}) > 8 $ indicate decisive evidence favoring the presence of a phase transition signal. As a complementary approach, the DIC incorporates both goodness-of-fit and model complexity through

where $ \bar{\theta} $ represents the posterior mean, $ D(\theta) = -2\ln{\cal{L}}(\theta) $, and $ p_D = \bar{D(\theta)} - D(\bar{\theta}) $ is the penalization term. The difference $ \Delta\text{DIC} = \text{DIC}_\text{null} - \text{DIC}_\text{FOPT} $ provides another measure of model preference, with larger positive values supporting the inclusion of the phase transition component.

Parameter estimation is performed using the nested sampling algorithm implemented in $ {\mathrm{dynesty}}$, accessed through the $ {\mathrm{Bilby}}$ Bayesian inference library. Figure 2 displays the resulting posterior distributions for a representative FOPT signal with amplitude $ \Omega_{\mathrm{PT}} = 3.9\times 10^{-11} $ and characteristic frequency $ f_{\mathrm{PT}} = 7\times 10^{-3} \mathrm{Hz} $. The recovered values, along with their median and $ 90 $% equal-tail uncertainties, are also summarized in Table 1.

Figure 2.  (color online) Posterior distributions of model parameters from Bayesian analysis using simulated Taiji data. The corner plot shows marginalized one-dimensional posteriors along the diagonal and joint two-dimensional distributions with $ 1\sigma $, $ 2\sigma $, and $ 3\sigma $ confidence contours in the off-diagonal panels. The red markers indicate the true parameter values used in generating the synthetic signal, which featured a FOPT with amplitude $ \Omega_{\mathrm{PT}} = 3.9\times 10^{-11} $ and characteristic frequency $ f_{\mathrm{PT}} = 7\times 10^{-3} \mathrm{Hz} $.

Our simulation framework incorporates a set of base parameters, including: the detector noise characterization parameters fixed at reference values of $ A = 3 $ and $ P = 8 $; Galactic foreground modeling with four parameters describing the DWD confusion noise: amplitude coefficients $ A_1 = 3.98 \times 10^{-16} $ and $ A_2 = 4.79 \times 10^{-7} $, with corresponding spectral slopes $ \alpha_1 = -5.7 $ and $ \alpha_2 = -6.2 $; ECB background parameterized by amplitude $ A_{\mathrm{ECB}} = 1.8 \times 10^{-9} $ with canonical spectral index $ \alpha_{\mathrm{ECB}} = 2/3 $. While these parameters remain constant throughout our analysis, it is important to note that each simulation represents a distinct statistical realization of the stochastic backgrounds, as the foreground components are characterized by their power spectral densities rather than deterministic waveforms.

Against this realistic background, we systematically inject phase transition signals spanning a two-dimensional parameter grid. The signal strength parameter $ \Omega_{\mathrm{PT}} $ and characteristic frequency $ f_{\mathrm{PT}} $ are varied across the following ranges:

This parameterization creates a grid of 100 distinct signal configurations, each requiring a separate Markov Chain Monte Carlo (MCMC) analysis. Figure 1 illustrates the frequency-domain representation of synthetic Taiji data for a representative case with $ \Omega_{\mathrm{PT}} = 3.9\times 10^{-11} $ and $ f_{\mathrm{PT}} = 7\times 10^{-3} \mathrm{Hz} $. The corresponding posterior distributions for this benchmark scenario are presented in Fig. 2, demonstrating that all model parameters are successfully recovered within the $ 2\sigma $ credible intervals.

Figure 3 and Fig. 4 display the measurement uncertainties in the recovered peak frequency $ f_{\mathrm{PT}} $ and amplitude $ \Omega_{\mathrm{PT}} $, respectively, across the parameter space. The error bars exhibit significant growth when $ \Omega_{\mathrm{PT}} $ falls below $ 1.4\times 10^{-11} $ or when $ f_{\mathrm{PT}} $ is less than $ 1.2 \times 10^{-3} \mathrm{Hz} $. This degradation in parameter estimation precision can be attributed to the competing influence of the DWD confusion background, which dominates the low-frequency sensitivity band of the detector and effectively masks cosmological signals below certain amplitude thresholds in this frequency regime.

Figure 3.  (color online) Comparison between injected and recovered peak frequencies ($f_\mathrm{PT}$) of the FOPT signal. Each point represents the median of the posterior distribution, with error bars indicating the 90% credible intervals. The dashed line represents perfect recovery.

Figure 4.  (color online) Comparison between injected and recovered amplitudes ($\Omega_\mathrm{PT}$) of the FOPT signal. Each point represents the median of the posterior distribution, with error bars indicating the 90% credible intervals. The dashed line represents perfect recovery.

Figure 5 presents the relative uncertainty in amplitude ($ \Delta\Omega_\mathrm{PT}/\Omega_\mathrm{PT} $), while Fig. 6 illustrates the relative uncertainty in peak frequency ($ \Delta f_\mathrm{PT}/f_\mathrm{PT} $) across the parameter space. As expected, $ \Delta\Omega_\mathrm{PT}/\Omega_\mathrm{PT} $ demonstrates a clear inverse relationship with signal strength, decreasing systematically as $ \Omega_\mathrm{PT} $ increases due to improved signal-to-noise ratio. Similarly, the fractional uncertainty in frequency determination $ \Delta f_\mathrm{PT}/f_\mathrm{PT} $ also diminishes with increasing signal amplitude. Notably, when the phase transition signal reaches $ \Omega_\mathrm{PT} \gtrsim 1.1\times 10^{-10} $, the frequency can be determined with high precision, achieving $ \Delta f_\mathrm{PT}/f_\mathrm{PT} \lesssim 0.1 $ across most of the frequency range.

Figure 5.  (color online) Measurement precision of the phase transition amplitude as a function of signal strength. The vertical axis shows the relative uncertainty ($\Delta\Omega_\mathrm{PT}/\Omega_\mathrm{PT}$) in the recovered amplitude, while the horizontal axis displays the injected amplitude values.

Figure 6.  (color online) Frequency resolution capabilities of the analysis pipeline across the detection band. The plot shows the relative uncertainty ($\Delta f_\mathrm{PT}/f_\mathrm{PT}$) in peak frequency estimation as a function of the injected signal frequency.

To quantitatively evaluate model selection capabilities, we present the logarithmic BFs in Fig. 7 and the DIC differences in Fig. 8, comparing models with and without the phase transition component across the parameter space. Both metrics exhibit consistent behavior, showing progressive improvement in detection confidence as $ \Omega_\mathrm{PT} $ increases. This concordance between independent statistical measures reinforces our confidence in the results. The observed trend aligns with theoretical expectations, as larger amplitude signals naturally produce more decisive evidence for the presence of a cosmological phase transition as opposed to the null hypothesis of only astrophysical and instrumental backgrounds.

Figure 7.  (color online) Model selection analysis using the Bayes factors. The heatmap displays logarithmic Bayes factors comparing models with and without a phase transition component, plotted as a function of signal amplitude ($\Omega_\mathrm{PT}$) and peak frequency ($f_\mathrm{PT}$).

Figure 8.  (color online) Model selection analysis using the Deviance Information Criterion (DIC). The heatmap illustrates the difference in DIC values between models with and without a phase transition component across the parameter space of signal amplitude ($\Omega_\mathrm{PT}$) and characteristic frequency ($f_\mathrm{PT}$).

Our comprehensive analysis demonstrated Taiji's significant potential for detecting and characterizing SGWBs from cosmological FOPTs. Through systematic Bayesian analysis incorporating realistic instrumental noise and astrophysical foregrounds, we found that Taiji can robustly detect phase transition signals with energy densities exceeding $ \Omega_{\mathrm{PT}} \gtrsim 1.4 \times 10^{-11} $ across most of its frequency band, with optimal sensitivity in the range of $ 10^{-3} $ to $ 10^{-2} $ Hz. For stronger signals with $ \Omega_{\mathrm{PT}} \gtrsim 1.1 \times 10^{-10} $, Taiji can determine the peak frequency with relative precision better than 10%. This sensitivity threshold represents a substantial improvement over current constraints [75, 76], enabling tests of various early universe scenarios, including strongly supercooled transitions and those associated with composite Higgs models or hidden sector physics [40, 77]. The consistency between our Bayesian evidence calculations and information-theoretic metrics provides a solid statistical foundation for future detection claims [78]. While our study focused on the broken power-law spectral template, future work should explore more physically motivated spectral shapes directly connected to specific phase transition parameters such as transition temperature, strength, and bubble wall velocity [22, 79].

Our analysis employed the standard approach of combining multiple TDI channels (A, E, T) from a single detector. While this method effectively suppresses instrumental noise through the null channel T, it exhibits inherent limitations for detecting stochastic backgrounds [80]. Single-detector analyses are fundamentally limited by the inability to distinguish between true GW signals and correlated instrumental artifacts. The null channel method, while useful for validation, cannot provide the same level of confidence as cross-correlation techniques between independent detectors.

For phase transition detection specifically, a multi-detector network could achieve detection thresholds potentially an order of magnitude lower than single-detector analyses, while providing more robust parameter estimation and reducing false positive rates. The different arm lengths and orientations of LISA (2.5 million km) and Taiji (3 million km) would offer complementary frequency responses, enhancing overall sensitivity across the millihertz band. The synergistic potential of Taiji operating concurrently with other space-based detectors, such as LISA, would further enhance detection prospects through cross-correlation techniques [81, 82].

We thank the anonymous referee for providing constructive comments and suggestions that greatly improve the quality of this manuscript. We acknowledge the use of HPC Cluster of ITP-CAS.