Symmetry Breaking in Parallel 1-K Sorption Coolers and Passive Suppression Strategy

1  Introduction

Sub-Kelvin cooling technology (T<1 K) has emerged as a fundamental prerequisite for frontier scientific exploration, playing an irreplaceable role in fields such as quantum computing and deep space exploration. A stable cryogenic environment is critical for suppressing thermal noise, thereby guaranteeing the precise operation of superconducting devices and high-sensitivity detectors. For instance, in Cosmic Microwave Background (CMB) observation missions—such as the Planck satellite or the upcoming LiteBIRD mission—bolometric detector arrays typically operate at 100 mK or lower to reach the photon shot noise limit [1–3]. Similarly, in superconducting quantum computing, the coherence time of qubits is directly constrained by environmental temperature, necessitating a rigorous operational environment of 10–20 mK [4–6].

As deep space exploration missions become increasingly complex, single-unit coolers struggle to meet escalating performance demands. To capture faint signals, next-generation astronomical telescopes plan to deploy Large-Scale Focal Plane Arrays (FPAs) containing thousands to tens of thousands of pixels, imposing stringent requirements on the cooling capacity at the 1-K stage [7–9]. Furthermore, specific long-duration observation missions require continuous temperature maintenance for days or even weeks. Although our previous work achieved extended hold times through optimized adsorption bed structures [10], the capacity limit of a single unit remains a bottleneck for overall system performance due to spatial constraints and efficiency limits. To overcome the capacity limitations of single units, a common engineering strategy is to operate multiple identical sorption coolers in parallel on a single thermal bus [11,12]. Theoretically, this modular parallel configuration should offer a linear superposition of cooling power: sharing the load between two coolers would halve the evaporation rate of each, significantly lowering the cold head temperature and doubling the hold time.

However, the stability of such parallel sub-Kelvin systems cannot be taken for granted. In traditional multiphase fluid networks, flow mal-distribution is a well-documented classical problem. For instance, based on the physical models established by Kakac and Bon [13], density-wave instabilities or dry-out phenomena in parallel channels are fundamentally hydrodynamic instabilities driven by the dynamic coupling of fluid flow resistance and latent heat. In contrast, sub-Kelvin parallel sorption coolers inherently lack macroscopic fluid flow between units; yet, in our dual-unit parallel experimental platform, we observed a previously unreported solid-state dynamic competition. Instead of synchronized cooling, the system exhibited “symmetry breaking”. The two coolers spontaneously evolved into an asymmetric mode where one unit captured the majority of the heat load and rapidly depleted its liquid helium, while the other remained dormant, causing the parallel system to degenerate into a “quasi-series” operation.

In-depth analysis reveals that the physical origin of this symmetry-breaking phenomenon lies in the strong non-linearity of thermal link materials and interfaces within the 0.8–1.5 K range. Governed by the Wiedemann-Franz law, the thermal conductivity of common thermal bus materials, such as Oxygen-Free High Thermal Conductivity (OFHC) copper, exhibits a linear temperature dependence (k(T)∝T) due to electron-impurity scattering [14]. Furthermore, based on Swartz’s Acoustic Mismatch Model (AMM) [15], the Kapitza thermal contact resistance (TCR) at solid-solid joints exhibits severe negative non-linearity (Rc∝T−3). This highly stochastic interface resistance often dominates the total thermal resistance below 4 K. When this material and interfacial non-linearity is coupled with the highly temperature-sensitive pumping dynamics of sorption coolers, it triggers a positive feedback loop: a marginal temperature elevation in one path sharply increases its thermal conductance, actively siphoning more heat flow from the common bus, which in turn raises the local temperature further, ultimately driving the system into a persistent thermal lock.

Existing research has primarily focused on the steady-state optimization of single units [16–18], leaving a critical research gap regarding transient thermodynamic instabilities in parallel sub-Kelvin configurations. To bridge this gap, this study systematically investigates the spontaneous symmetry-breaking mechanism and its suppression. The main contributions are:

(i)   Phenomenological Reproduction: We experimentally quantify the reproducible degeneration from “parallel” to “quasi-series” operation.

(ii)   Mechanism Elucidation: Distinct from classical fluid instabilities, we derive a theoretical stiffness-competition framework, mathematically explaining how solid-state heat conduction drives the “avalanche-like” heat flow bifurcation.

(iii)   Load Dependence & Passive Suppression: We reveal that higher thermal loads paradoxically intensify the destabilizing positive feedback. Consequently, we propose a ballast thermal resistance concept as a promising, numerically supported strategy to mitigate the asymmetry and offer valuable design guidelines for recovering parallel cooling performance.

2  Experimental Setup

2.1 Cryogenic Platform Description

The schematic of the parallel sorption cooling system is illustrated in Fig. 1. The experimental platform utilizes a closed-cycle Helium-4 sorption refrigeration cycle to achieve sub-Kelvin temperatures.

Figure 1: Schematic diagram of the parallel sorption cooling system.

The backbone of the system is a two-stage G-M cryocooler (Model: Sumitomo RDK-101D cold head coupled with an HC-4F compressor unit), which provides the necessary pre-cooling power. This cryocooler delivers cooling capacities of approximately 0.1 W at 4.2 K and 5 W at 40 K, serving as the heat sink for the sorption coolers. The entire assembly is housed within a custom-manufactured stainless steel vacuum dewar. To minimize convective heat transfer, a molecular pump unit maintains a high-vacuum environment with a background pressure of approximately 1.0×10−5 Pa.

The physical assembly of the system is depicted in Fig. 2. Two technically identical sorption coolers (labeled Unit A and Unit B) are installed in parallel on the 4 K cold plate, as shown in Fig. 2a. Considering the demanding requirements for vacuum integrity and structural stability, all critical joints between copper and stainless steel components were realized through vacuum brazing or TIG welding.

Figure 2: Photographs of the experimental setup: (a) two identical helium sorption coolers (HSC) mounted on the Sumitomo G-M cryocooler; (b) detailed view of the parallel OFHC copper thermal connection between the two cold heads and the shared busbar.

The evaporator cold heads of the two units are thermally linked to a shared 1-K busbar (Payload Platform) via Oxygen-Free High Thermal Conductivity (OFHC) copper straps. The specific connection of this thermal link, which is identified as the source of the symmetry breaking, is shown in Fig. 2b.

To suppress radiative heat leaks, a two-stage radiation shield system is integrated, anchored to the first (40 K) and second (4 K) stages of the G-M cryocooler, respectively. Both shields are wrapped with 30 layers of Multi-Layer Insulation (MLI), composed of 20 μm thick aluminized polyimide film with polyester spacers. Furthermore, to minimize conductive heat leaks, the 1-K platform is mechanically suspended from the 4-K stage using Kevlar ropes, effectively isolating it from vibration and thermal conduction.

2.2 Measurement and Data Acquisition

The system employs a rigorous thermometry setup to capture transient dynamics. Cernox resistance temperature sensors (Lake Shore Cryotronics, Model CX-1050) are mounted on the 1-K platform and the cold heads of both units (indicated by blue squares in Fig. 1). These sensors are calibrated for the 0.3–300 K range, offering high sensitivity in the sub-Kelvin region.

The thermal load on the 1-K platform is simulated using a ceramic resistive heater. The heating power (Qload) is applied via a precision DC power source and measured using the four-wire method to eliminate lead resistance errors.

All temperature and power signals are acquired by a Lake Shore 350 temperature controller and recorded by a host computer at a sampling frequency of 1 Hz. This rate is sufficient to capture the thermal relaxation processes during the startup and locking phases.

2.3 Experimental Protocol and Uncertainty

To investigate the symmetry breaking phenomenon, a variable load experiment was designed. The heating power applied to the 1-K platform was gradually increased from 0 to 5 mW.

Definition of Start Time (ts,i): Defined as the moment when the temperature curve of the cold head exhibits a sharp inflection point, accompanied by a significant increase in the cooling rate of the platform.

Definition of Depletion (tdepletion): To ensure the strict reproducibility of our analysis, the end of the “effective cooling state” (i.e., the moment of macroscopic liquid helium depletion) is defined using a mathematical threshold. It is recorded as the exact timestamp when the time derivative of the cold head temperature (dTc/dt) consecutively exceeds 10 mK/s. It is important to note that the extracted depletion time is highly insensitive to the specific numerical choice of this threshold. Prior to depletion, the steady-state temperature drift is minimal (e.g., typically <24 mK over a 30-min period). Upon liquid helium exhaustion, the loss of latent heat cooling triggers a severe step-like phase change, causing the local temperature to surge by approximately 240 mK within 10 s. This stark physical contrast ensures that the depletion point is unambiguously identifiable.

Uncertainty Analysis: The experimental uncertainty primarily stems from sensor calibration and instrument precision. The specifications of the Cernox sensors are listed in Table 1, demonstrating a precision of better than ±4 mK within the primary operational range of 0.8–1.5 K. It is critical to note that the observed macroscopic asymmetry significantly overshadows this measurement error. For instance, under a typical 2 mW load, the temperature jump upon helium depletion reaches approximately 400 mK, and the steady-state temperature differential between the active and depleted evaporators exceeds 180 mK. Because this physically driven temperature deviation is nearly two orders of magnitude larger than the sensor precision, the propagation of sensor uncertainty has a negligible impact on the computed Synchronization Ratio and the fundamental conclusions regarding thermal stability. The combined relative uncertainty of the applied heat load is estimated to be within ±0.5%.

Experimental Duration and Statistical Reliability: Given the complex nature of sub-Kelvin cryogenics, a single complete operational cycle of the parallel sorption cooler—encompassing physical adsorption, Helium-4 condensation, adiabatic pumping, and the subsequent thermal recovery phase—requires approximately 4 days. To ensure the statistical reliability of the observed thermodynamic instabilities, 4 completely independent repeated experimental runs were conducted under identical macroscopic configurations, representing an accumulated ultra-low temperature operational time of over 16 days. Across all independent trials, the symmetry-breaking phenomenon exhibited a 100% occurrence rate. While the precise transient temperature trajectories displayed marginal variations between runs—attributed to the extreme sensitivity to initial conditions inherent in strongly non-linear thermal networks—the fundamental physical trend remained consistently reproducible: the system spontaneously degenerated into an asymmetric mode.

Furthermore, for the load-dependence analysis, experiments were conducted across six discrete, continuously increasing active loads (0, 1, 2, 3, 4, 5 mW). The system exhibited a strictly monotonic degradation trend in synchronization as the load increased. From a statistical perspective, the emergence of a monotonic macroscopic trajectory across six discrete states indicates that the observed degeneration is driven by a deterministic thermodynamic mechanism, rather than random experimental dispersion.

3  Theoretical Framework and Numerical Model

To quantitatively investigate the symmetry breaking mechanism and evaluate the proposed suppression strategy, a transient thermal-fluidic resistance network model was developed. This model abstracts the complex physical system into a network of discrete thermal nodes connected by non-linear thermal resistors, as schematically shown in Fig. 3.

Figure 3: Schematic diagram of the equivalent circuit for the thermal resistance network.

3.1 Thermal-Fluidic Resistance Network

To simplify the complex physical system, the model is discretized into three primary nodal categories: the payload platform (Tp), the two parallel cold heads (Tc,1, Tc,2), and the environmental heat sink (Tenv). The formulation is based on the Lumped Capacitance Method.

This lumped-parameter assumption is rigorously justified by evaluating the Biot number (Bi=hLc/k). The liquid helium pool and the high-purity OFHC copper structures within the cold head exhibit thermal conductivities that are sufficiently high to ensure that the internal conduction resistance is mathematically negligible compared to the thermal resistance of the external connecting copper straps and the solid-solid interfaces. A quantitative calculation for the discrete nodes confirms that Bi≪0.1, validating the assumption of uniform temperature distribution within each component. Furthermore, the activated carbon inside the sorption pump is tightly compacted, yielding a minimal characteristic conduction length. Consequently, internal temperature gradients within the active components are treated as mathematically negligible during transient heat transfer processes.

A critical feature of this model is the incorporation of temperature-dependent thermophysical properties. The thermal conductivity of the OFHC copper linkages (kCu) is not constant but varies significantly with temperature.

Although represented as a single equivalent resistor RCu in the schematic (Fig. 3), in the numerical solver, this resistance is computed iteratively by discretizing the link into n segments. This integration accounts for the temperature gradient along the strap based on the geometric parameters and the material properties:

RCu(T)=∫0LdxA⋅k(T(x)),(1)

where A and L denote the cross-sectional area and effective length of the copper strap, respectively. The thermal conductivity k(T) is derived from the NIST cryogenic material database. This non-linear term, which increases with temperature in the 1–10 K range, is identified as the source of the positive feedback loop.

The transient thermal evolution is governed by the energy conservation principle. For a generic node i with total heat capacity Ci (including addenda), the instantaneous energy balance is described by the Ordinary Differential Equation (ODE):

CidTidt=∑jGij⋅(Tj−Ti)+Qin,i,(2)

where Ti and Tj are the instantaneous temperatures of the node and its neighbor, Gij is the thermal conductance between them, and Qin,i represents external heat sources.

Applying this to the 1-K Platform (Tp):

CpdTpdt=−∑i=1,2Gi(Tavg)⋅(Tp−Tc,i)−Gl(Tp−Tenv)+Qload,(3)

where Cp is the platform heat capacity, Gi is the conductance of the link to the i-th cold head, Gl represents the parasitic heat leak (radiation and support conduction), Tenv is the fixed environmental temperature (4.2 K), and Qload is the experimental heating power.

For the Cold Head nodes (Tc,i):

Cc,idTc,idt=Gi(Tavg)⋅(Tp−Tc,i)−Qcool,i(Tc,i),(4)

where Cc,i is the lumped heat capacity of the cold head and liquid helium pool. The term Qcool,i represents the net cooling power provided by evaporation, modeled as a function of the saturation vapor pressure and the latent heat of Helium-4, constrained by the finite liquid inventory.

3.2 Theoretical Stability Criterion via Stiffness Competition

To theoretically define the stability boundaries of the system, this section investigates the robustness of equilibrium points by analyzing the steady-state energy balance of a single branch using perturbation analysis.

The operational state of a cold head (Tc) is determined by the intersection of the Heat Supply Curve (Qin) and the Cooling Capacity Curve (Qcool).

Heat flows from the common platform (Tp) to the individual cold head (Tc). In the 0.8–1.5 K regime, phonon-mediated heat conduction (∝T3) is effectively frozen out; therefore, heat transfer in high-purity OFHC copper is primarily governed by electron-impurity scattering [19]. Governed by the Wiedemann-Franz law, the thermal conductivity exhibits a strong linear temperature dependence (k(T)≈βT). Furthermore, the thermal resistance of these external OFHC links and solid interfaces is orders of magnitude larger than the internal evaporation thermal resistance of the liquid helium pool. Thus, the non-linearity of the external conduction network fundamentally constitutes the primary thermal bottleneck and dictates the heat flow distribution. Integrating Fourier’s law yields:

Qin(Tc)=AL∫TcTpβTdT=Aβ2L(Tp2−Tc2),(5)

where A/L is the geometric factor and β is the temperature coefficient of conductivity. It is important to emphasize that β is not an arbitrary mathematical fitting parameter; rather, it is a physical constant derived directly from the NIST cryogenic material database. Curve Topology: This equation describes a downward-opening parabola. As Tc and Tp rise, even if the driving temperature difference remains constant, the thermal conductivity increases significantly. This results in a non-linear surge in heat flow within the specific temperature range.

The cooling power is governed by the mass flow rate driven by the saturation vapor pressure. Combining the flow equation with the Clausius-Clapeyron relation, we obtain:

Qcool(Tc)=η⋅[P0exp⁡(−BTc)−Ppump],(6)

where η is a coefficient related to the flow conductance and latent heat, and B is the specific gas constant for Helium-4. Eq. (6) is derived from classical evaporative cooling gas dynamics. The parameter η aggregates the effective pumping speed and the latent heat of evaporation. These specific parameters were strictly determined through independent single-unit steady-state characterization tests prior to the parallel experiments. This represents a convex exponential growth curve. It remains flat in the low-temperature region and rises steeply as temperature increases.

A mathematical intersection of these curves (Qin=Qcool) is not necessarily a physically stable operating point. For a point to be stable, the system must possess a net restoring force. This requires that the “stiffness” of the cooling mechanism exceeds the “gain” of the thermal link. We define two stiffness parameters:

Destabilizing Gain (Kcond): The magnitude of the slope of the heat input curve, representing the positive feedback strength of the thermal link.

Kcond=|dQindTc|=AβLTc,(7)

Restoring Stiffness (Kcool): The slope of the cooling curve, representing the negative feedback capability.

Kcool=dQcooldTc≈ηP0BTc2exp⁡(−BTc),(8)

Based on the criterion Kcond>Kcool (instability condition), we can explain the observed behaviors through scaling analysis:

Operational Window (Instability Zone): In the experimental temperature range (approx. 0.8–1.5 K), the thermal conductivity of copper rises sharply (Kcond∝Tc). Our analysis suggests that in this specific window, the increase in Destabilizing Gain outpaces the increase in Restoring Stiffness. Consequently, higher loads (which push Tc higher) intensify the positive feedback. If a minute temperature perturbation occurs, the resulting heat-siphoning effect overwhelmingly exceeds the cooler’s responsive capacity, theoretically triggering the observed “avalanche-like” heat flow bifurcation and symmetry breaking.

High-Temperature Limit (Restoration): It is important to note that the positive feedback does not increase infinitely. At sufficiently high temperatures (e.g., >2 K), the exponential nature of the cooling power term (Kcool) will eventually overtake the linear growth of the thermal conductivity (Kcond). However, within the nominal working range of the sorption cooler, the material non-linearity dominates, leading to the observed symmetry breaking.

3.3 Principle of Ballast Thermal Resistance

According to the derivation in Section 3.2, the root cause of thermal locking is an excessive Destabilizing Gain in the thermal link (Kcond>Kcool). To restore stability without altering the cooler’s intrinsic physics (i.e., maintaining Kcool), the engineering solution is to suppress Kcond. We propose inserting a temperature-insensitive Ballast Thermal Resistance (Rb), such as a stainless steel or brass shim, in series between the cold head and the payload.

With the ballast resistor, the total thermal resistance Rtot becomes the sum of the non-linear copper resistance RCu(Tc) and the constant resistance Rb:

Rtot(Tc)=RCu(Tc)+Rb,(9)

The modified heat flow Qin′ is now expressed as:

Qin′(Tc)=Tp−TcRCu(Tc)+Rb,(10)

To rigorously evaluate the system’s stability, we differentiate Eq. (10) with respect to Tc using the quotient rule. Recognizing that ∂(Tp−Tc)/∂Tc=−1 and ∂Rb/∂Tc=0, the derivative yields:

∂Qin′∂Tc=−1RCu+Rb−Tp−Tc(RCu+Rb)2∂RCu∂Tc,(11)

This mathematically corrected derivative reveals a profound physical mechanism, decomposing the thermal response into two competing effects:

1.   Linear Restoring Term (First Term): This term is negative (−1RCu+Rb). Physically, it indicates that an elevation in cold head temperature naturally decreases the driving temperature difference (Tp−Tc), thereby reducing the incoming heat flow. This serves as an inherent negative feedback (stabilizing) mechanism.

2.   Non-Linear Destabilizing Term (Second Term): Because the absolute thermal resistance of OFHC copper drops sharply as temperature rises in this regime (∂RCu/∂Tc<0), this entire second term becomes positive. This represents the “black hole” effect—the avalanche-like positive feedback where a local temperature rise drastically reduces resistance and siphons additional heat flow.

The modified Destabilizing Gain (Kcond′) is defined by the magnitude of this positive feedback term:

Kcond′=|Tp−Tc(RCu+Rb)2∂RCu∂Tc|,(12)

Comparing this to the original system (Rb=0), the introduction of Rb exerts a powerful, multi-layered stabilizing effect:

First, it dilutes the overall weight of the copper’s thermal fluctuation. More importantly, it squares the denominator of the destabilizing term—transforming it to (RCu+Rb)2—thereby suppressing the non-linear positive feedback at a quadratic rate. By selecting a sufficiently large Rb, the constant resistance dictates the overall heat transfer, allowing the linear restoring behavior to overwhelm the non-linear avalanche effect, effectively enforcing the condition Kcond′<Kcool and restoring uniform load distribution.

3.4 Numerical Implementation and Initial Conditions

Due to the dramatic variation in the thermal conductivity of copper at low temperatures, the governing equations exhibit highly non-linear characteristics. The simulation was implemented in MATLAB using the ode45 solver (based on an explicit Runge-Kutta (4,5) formula) with an adaptive step-size algorithm. The relative tolerance and absolute tolerance were set to 1×10−6 and 1×10−9, respectively, to accurately capture the rapid thermal transient processes.

To reproduce the experimentally observed symmetry breaking, the mathematical symmetry of the ideal model must be broken. At t=0, a minute initial geometric asymmetry was introduced as a perturbation source. Specifically, the geometric factor (A/L) of one thermal branch was set to have a 0.5% deviation (A1=0.995A2). This specific value is physically grounded: a 0.5% deviation directly corresponds to a standard precision machining tolerance of ±50µm on a typical 10-mm OFHC copper linkage. Beyond macroscopic manufacturing limits, this deterministic parameter serves as an “equivalent lumped perturbation” in our network model. It mathematically represents the inherent, unmeasurable microscopic stochasticity in real cryogenic systems, such as the 10%–20% natural fluctuation in Kapitza Thermal Contact Resistance (TCR) across mechanically identical joints, as well as the minute gas-dynamics delays during the simultaneous startup of the sorption pumps. This setup allows the numerical model to deterministically capture the system’s sensitive dependence on initial conditions and reliably predict the bifurcation behavior without over-fitting unmeasurable random variables.

To ensure the rigor of the numerical validation, it is necessary to explicitly distinguish between independently measured physical constants and assumed modeling parameters. Table 2 summarizes the classification and origin of the key parameters utilized in the transient thermal-fluidic model. Notably, the parameters dictating the intrinsic cooling capacity (e.g., Tbase, η, B) were determined via independent single-unit characterization tests, prior to the assembly of the parallel network. The model does not rely on post-hoc curve fitting of the parallel transient data; rather, the macroscopic bifurcation behavior is predictively driven by the established material non-linearities and the equivalent lumped perturbation.

4  Results and Discussion

4.1 Phenomenological Asymmetry

To characterize the transient operational dynamics of the parallel sorption cooling system, cooldown experiments were conducted under varying thermal loads. Fig. 4 presents the temporal evolution of the two cold heads (Tc,1,Tc,2) under a zero-load condition (0 mW), where the system is subject only to parasitic heat leaks from the environment.

Figure 4: Temporal evolution of cold head temperatures under 0 mW load (parasitic load only).

Initially, the two coolers exhibit a high degree of synchronization during the rapid cooldown phase. However, as the temperature drops below 0.9 K, a symmetry breaking becomes observable: Tc,1 stabilizes at a slightly higher value than Tc,2. This reveals a critical physical reality: despite the symmetric design, parasitic heat leaks (from 4 K radiation and sensor wiring) are not perfectly distributed, imposing a marginally higher load on Unit 1. Although subtle, this imbalance accumulates over time. Unit 1 depletes its liquid helium inventory first at t≈12 h, causing its temperature to rise sharply. Consequently, the remaining parasitic load is transferred entirely to Unit 2, resulting in a distinctive “step-like” temperature increase for the second unit. Despite this sequential depletion, both units operate simultaneously for the majority of the hold time. Statistical analysis indicates that the “relay” phase (where Unit 2 operates alone) accounts for only 23% of the total duration, suggesting that the system maintains a relatively high degree of parallelism under low-load conditions.

Subsequent experiments were conducted under an active thermal load of 3 mW, yielding significantly altered thermodynamic behavior, as shown in Fig. 5. First, the equilibrium temperature of the cold heads is noticeably higher than in the 0 mW case, pushing the OFHC copper linkages into a temperature regime with significantly higher thermal conductivity. While both units share the load during the initial phase, the asymmetry in heat flow distribution is drastically amplified. Unit 1 rapidly dominates the heat absorption, depleting its liquid helium much faster than Unit 2.

Figure 5: Temporal evolution of cold head temperatures (solid lines) and adsorption pump temperatures (dashed lines) under a 3 mW active load.

This phenomenon is strongly corroborated by the temperature profiles of the adsorption pumps (indicated by dashed lines in Fig. 5). The temperature rise of the adsorption pump serves as a direct proxy for the adsorption rate (and thus the instantaneous cooling power), as the adsorption process is exothermic. During the first 2 h, the temperature of Pump 1 is significantly higher than that of Pump 2, confirming that Unit 1 is handling the majority of the mass flow and cooling load. Conversely, once Unit 1 is depleted (indicated by the drop in Pump 1’s temperature), Pump 2’s temperature rises sharply as it takes over the load. This confirms that under higher loads, the system degenerates from a “parallel” mode into a pronounced “relay” or “quasi-series” mode.

To quantify the degree of parallelism across different operating conditions, we define the Synchronization Ratio (S) as:

S=toverlapttotal,(13)

where toverlap represents the duration during which both coolers are in an effective cooling state (i.e., temperatures remaining below the strictly defined dT/dt>10 mK/s depletion threshold), and ttotal is the total effective hold time of the 1-K platform. A value of S close to 1 indicates ideal parallel operation, while a lower value implies a tendency towards serial or relay operation.

Fig. 6 illustrates the variation of the Synchronization Ratio S with respect to the applied load power. The experimental data reveals a significant physical trend: the synchronization ratio decreases markedly as the thermal load increases. Low Load (0 mW): The system achieves S=75.3, indicating that both units contribute to cooling for the majority of the time. High Load (3 mW): The ratio drops to 51.3, implying that for nearly half the operation time, the system is effectively running on a single unit.

Figure 6: Total hold time and synchronization ratio (S) as a function of applied heating power.

This inverse relationship aligns with the theoretical “Stiffness Competition” mechanism proposed in Section 3. In high-load (high-temperature) scenarios, although the restoring stiffness of the cooler (Kcool) increases with vapor pressure, the destabilizing gain of the OFHC copper link (Kcond) increases more aggressively due to material non-linearity. The reduced thermal resistance of the copper straps at higher temperatures allows minute initial temperature differences to drive massive disparities in heat flow distribution. In other words, at high loads, the material-induced destabilizing gain overwhelms the restoring stiffness, locking the system into an asymmetric state. Conversely, at low loads, the copper straps operate in a lower conductivity regime, reducing the positive feedback strength and allowing the system to maintain synchronization more easily.

4.2 Numerical Model Verification

To validate the theoretical mechanism proposed in Section 3, we reproduced the experimental scenarios using the transient numerical model. Considering the inevitable manufacturing tolerances and variations in contact thermal resistance inherent in the experimental setup, a geometric asymmetry of 0.5% was introduced into the model as an initial condition. Fig. 7 illustrates the simulated transient evolution of cold head temperatures and heat flow distribution under a 1000 µW load. The simulation clearly reconstructs the “thermal locking” process driven by positive feedback. At the instant of startup, the heat flow is distributed almost equally between the two branches. However, as the platform temperature rises slightly, the branch with marginally lower thermal resistance (Unit 1) experiences a faster increase in thermal conductivity. This reduces its resistance further, attracting a larger fraction of the heat flow in an “avalanche-like” process that rapidly drives the system into a locked state where Unit 1 shoulders nearly 100% of the load. By overlaying the experimental temperature data as scatter points in Fig. 7a, we observe that the model exhibits excellent agreement with the measurements, accurately capturing both the bifurcation timing and the magnitude of the steady-state temperature differences. To further ensure the rigorousness of this validation, quantitative statistical indicators were introduced. Recognizing that highly non-linear avalanche processes are extremely sensitive to time-domain shifts, we evaluated the macroscopic thermodynamic features rather than relying solely on point-by-point time-series metrics. The steady-state base temperature predicted by the model exhibits an absolute error of approximately 10 mK (a relative error of merely ~1%) compared to the experimental measurements. Furthermore, the Mean Absolute Error (MAE) of the Synchronization Ratio (S) across different loads is less than 0.06. These metrics robustly support the model’s fidelity in capturing both the bifurcation timing and the magnitude of the steady-state load imbalance.

Figure 7: Transient response and verification under a 1000 µW load. (a) Comparison of simulated and experimental cold head temperatures. (b) Simulated heat flow distribution between the two parallel units over time.

Beyond the transient behavior, we further validated the model’s accuracy across the full operational spectrum by examining the steady-state thermal response. Fig. 8 compares the experimental and simulated steady-state temperatures of the cold heads under varying thermal loads (0–5 mW). The deviation between the model predictions and experimental measurements remains consistently below 20 mK. Crucially, the cold head temperature exhibits a quasi-linear increase with the applied load. While this trend appears simple visually, it implies a highly non-linear evolution in load distribution. As established in Section 3, the cooling capacity (Qcool) increases exponentially with temperature. Therefore, a linear rise in steady-state temperature indicates that the actual cooling power provided by the dominant cold head is increasing exponentially. This confirms that as the total load increases, the dominant head assumes a progressively larger proportion of the heat load (i.e., the load share increases non-linearly), thereby exacerbating the system’s imbalance.

Figure 8: Steady-state temperature validation.

Based on the validated model, we investigated the sensitivity of the system to manufacturing precision to determine if the symmetry breaking could be mitigated by tighter tolerances. Fig. 9 displays the simulated heat flow evolution under varying degrees of initial geometric asymmetry (ϵ), ranging from 0.01% to 0.5%. The results indicate that the system operates in a state of critical instability. Reducing the asymmetry from 0.5% to 0.1%—a level far exceeding typical machining tolerances—only delays the onset of the positive feedback by approximately 10 s. Once the bifurcation begins, the positive feedback loop inevitably drives the system to the same asymmetric steady state, with a temperature difference of less than 1 mK between the cases. This finding theoretically rules out improving manufacturing precision as a viable engineering solution, thereby reinforcing the necessity for the passive ballast resistance strategy proposed in the following section.

Figure 9: Sensitivity analysis: transient evolution of heat flow distribution under different initial geometric asymmetries (ϵ).

4.3 Optimization of Passive Suppression Strategy

To resolve the thermal locking issue without increasing system complexity, we propose a passive suppression strategy based on a “Ballast Thermal Resistance” (Rb). The core principle involves inserting a temperature-insensitive constant thermal resistor (e.g., a stainless steel or brass shim) in series between the cold head and the load. This linear resistance serves to dilute and suppress the positive feedback gain induced by the non-linearity of the copper straps.

To quantify the efficacy of this strategy and identify the optimal resistance value, Fig. 10 illustrates the variation of system performance parameters with respect to Rb under a 3 mW load. The plot utilizes a dual-axis format to juxtapose the benefit—Synchronization Ratio (S) on the left axis—against the cost—Platform Temperature Penalty (ΔTpenalty) on the right axis.

Figure 10: Trade-off analysis: synchronization ratio (S) and platform temperature penalty (ΔTpenalty) as a function of ballast thermal resistance (Rb).

Analyzing the recovery of synchronization, the system undergoes a distinct stability phase transition as Rb increases from zero:

Strong Locking Zone (Rb<2 K/W): In this region, the additional damping is insufficient to counteract the strong positive feedback gain of the copper straps at high temperatures, keeping the synchronization ratio at a low level.

Transition Zone (2–10 K/W): As Rb increases, the system enters a rapid recovery phase. The synchronization ratio S climbs sharply in an S-shaped curve, indicating that the heat flow distribution between the two branches is quickly converging towards equilibrium.

Saturation Zone (Rb>10 K/W): Once Rb exceeds a critical threshold, the synchronization index saturates and stabilizes above 0.95, signifying that the asymmetric bifurcation has been effectively suppressed and synchronous startup characteristics are fully restored.

However, introducing extra thermal resistance inevitably increases the heat transfer temperature difference, raising the minimum achievable temperature of the 1-K platform. The right-axis curve shows that the temperature penalty (ΔTpenalty) increases quasi-linearly with Rb. By weighing the stability benefits against the thermal penalty, a clear law of diminishing returns is observed: beyond 10 K/W, further increases in resistance yield negligible improvements in synchronization (S is already approaching 1.0), while the temperature penalty continues to accumulate.

Consequently, we identify Rb≈10 K/W as the optimal engineering design point for this system. At this resistance value, the system achieves a synchronization probability exceeding 95% with a minimal temperature penalty of less than 20 mK (raising the platform temperature from 1.00 to 1.02 K). This demonstrates that the proposed passive strategy can effectively enforce uniform load distribution between the two parallel HSCs with minimal performance sacrifice.

Practical Implementation and Reliability: From a practical engineering perspective, translating this theoretical ballast resistance into a physical assembly requires careful material selection. We recommend the insertion of commercial Stainless Steel (e.g., 304/316) or Brass shims. In the 1-K regime, the thermal conductivity of stainless steel is orders of magnitude lower than that of OFHC copper and exhibits a significantly weaker temperature dependence. This characteristic provides a relatively constant, predictable, and “linear” thermal resistance that effectively dampens the non-linear positive feedback of the copper straps. Furthermore, commercial shims ensure excellent manufacturability with precise thickness control at minimal cost. However, mechanical reliability must be carefully managed. The inherent mismatch in the Coefficient of Thermal Expansion (CTE) between the copper strap and the stainless steel shim from 300 K down to 1 K can cause standard bolted joints to loosen, leading to an unpredictable surge in Thermal Contact Resistance (TCR). To ensure thermomechanical stability, it is highly recommended to integrate Beryllium-Copper (Be-Cu) Belleville washers into the joint assembly. These spring washers provide a constant mechanical preload that compensates for thermal contraction, ensuring reliable interfacial contact.

Scalability to Multi-Cooler Arrays (N>2): When extending this stabilization strategy to larger systems with N parallel coolers, the stability requirements become increasingly stringent. If a temperature fluctuation occurs in one unit, it faces the combined heat-siphoning capacity of the remaining (N−1) active coolers acting as a massive heat sink, thereby intensifying the non-linear avalanche effect. For such multi-cooler configurations, we suggest arranging the ballast resistance in a “Star-network topology”. It is structurally insufficient to place a single resistor on the main thermal bus; instead, an independent, matched ballast resistor (shim) should be connected in series within every individual branch linking a cold head to the common payload. This topological design establishes physical thermal isolation between parallel branches, mitigating the cross-talk positive feedback loop.

5  Conclusions

This study systematically investigated the transient thermal stability of parallel Helium-4 sorption coolers, revealing a “symmetry breaking” phenomenon where the system spontaneously degenerates into an asymmetric quasi-series mode. Through theoretical modeling and experimental validation, we elucidated that this thermal locking is driven by a positive feedback loop originating from the temperature-dependent thermal conductivity (k∝T) of the OFHC copper linkages, which momentarily overwhelms the restoring stiffness of the sorption cooling mechanism.

Based on these findings, we outline the following design guidelines for parallel sub-Kelvin architectures:

1.   Failure of Uniform-Flow Assumptions: In the sub-Kelvin regime, traditional macroscopic uniform-flow assumptions may be invalid. The strong non-linearity of solid-state conduction materials and contact resistances should be integrated into the dynamic coupling analysis.

2.   Stiffness Criterion: The stability of the parallel network is governed by the mathematical competition between the heat-siphoning “Destabilizing Gain” of the thermal linkages and the “Restoring Stiffness” of the sorption coolers.

3.   Passive Suppression via Star Topology: Merely improving manufacturing precision is insufficient to prevent instability. We suggest incorporating low-thermal-conductivity materials (such as stainless steel shims equipped with Belleville washers) as ballast resistors. For multi-unit systems (N>2), employing a star-network topology—where independent ballast resistors are installed on each branch—is a promising, numerically supported strategy to enforce uniform load distribution at the cost of an acceptable absolute temperature penalty.

Finally, we acknowledge certain limitations of the current study. The stability criterion and numerical model primarily address quasi-steady thermal loads and do not fully capture the transient response under extreme pulsed heat loads (e.g., instantaneous detector calibration pulses). Additionally, the lumped-parameter approach simplifies the complex, non-equilibrium gas-dynamics and mass transfer within the porous activated carbon bed. Addressing these dynamic complexities presents a valuable direction for future multi-dimensional cryogenic thermal-fluidic modeling.

Acknowledgement: Not applicable.

Funding Statement: This work is supported by the National Natural Science Foundation Projects (52576028), the Hundred Talents Program of the Chinese Academy of Sciences, the Strategic Priority Research Program of Chinese Academy of Sciences (XDB35000000, XDB35040102).

Author Contributions: The authors confirm contribution to the paper as follows: conceptualization, Shaoshuai Liu and Yinong Wu; methodology, Lihao Lu; software, Lihao Lu; validation, Yan Lu, Zhenhua Jiang and Yinong Wu; formal analysis, Lihao Lu; investigation, Lihao Lu; resources, Shaoshuai Liu; data curation, Lihao Lu; writing—original draft preparation, Lihao Lu; writing—review and editing, Shaoshuai Liu and Yan Lu; visualization, Lihao Lu, Yan Lu and Zhenhua Jiang; supervision, Shaoshuai Liu; project administration, Shaoshuai Liu; funding acquisition, Shaoshuai Liu. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the Corresponding Author, Shaoshuai Liu, upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Hazumi M, Ade PA, Adler A, Allys E, Arnold K, Auguste D, et al. LiteBIRD satellite: JAXA’s new strategic L-class mission for all-sky surveys of cosmic microwave background polarization. In: Lystrup M, Batalha N, Tong EC, Siegler N, Perrin MD, editors. Space telescopes and instrumentation 2020: optical, infrared, and millimeter wave. Bellingham, WA, USA: SPIE; 2020. [Google Scholar]

2. Bhandari P, Prina M, Bowman RC, Paine C, Pearson D, Nash A. Sorption coolers using a continuous cycle to produce 20 K for the Planck flight mission. Cryogenics. 2004;44(6–8):395–401. doi:10.1016/j.cryogenics.2004.02.022. [Google Scholar] [CrossRef]

3. Gualtieri R, Filippini JP, Ade PAR, Amiri M, Benton SJ, Bergman AS, et al. SPIDER: CMB polarimetry from the edge of space. J Low Temp Phys. 2018;193(5):1112–21. doi:10.1007/s10909-018-2078-x. [Google Scholar] [CrossRef]

4. Arute F, Arya K, Babbush R, Bacon D, Bardin JC, Barends R, et al. Quantum supremacy using a programmable superconducting processor. Nature. 2019;574(7779):505–10. doi:10.1038/s41586-019-1666-5. [Google Scholar] [PubMed] [CrossRef]

5. Krantz P, Kjaergaard M, Yan F, Orlando TP, Gustavsson S, Oliver WD. A quantum engineer’s guide to superconducting qubits. Appl Phys Rev. 2019;6(2):021318. doi:10.1063/1.5089550. [Google Scholar] [CrossRef]

6. Blais A, Grimsmo AL, Girvin SM, Wallraff A. Circuit quantum electrodynamics. Rev Mod Phys. 2021;93(2):025005. doi:10.1103/revmodphys.93.025005. [Google Scholar] [PubMed] [CrossRef]

7. Galitzki N, Ali A, Arnold KS, Ashton PC, Austermann JE, Baccigalupi C, et al. The simons observatory: instrument overview. In: Zmuidzinas J, Gao JR, editors. Millimeter, submillimeter, and far-infrared detectors and instrumentation for astronomy IX. Bellingham, WA, USA: SPIE; 2018. [Google Scholar]

8. Holland WS, Bintley D, Chapin EL, Chrysostomou A, Davis GR, Dempsey JT, et al. SCUBA-2: the 10,000 pixel bolometer camera on the James Clerk maxwell telescope. Mon Not R Astron Soc. 2013;430(4):2513–33. doi:10.1093/mnras/sts612. [Google Scholar] [CrossRef]

9. Sobrin JA, Ade P, Ahmed Z, Anderson A, Avva J, Basu Thakur R, et al. Design and characterization of the SPT-3G receiver. In: Zmuidzinas J, Gao JR, editors. Millimeter, submillimeter, and far-infrared detectors and instrumentation for astronomy IX. Bellingham, WA, USA: SPIE; 2018. [Google Scholar]

10. Liu T, Cui X, Lu L, Ying K, Wang Y, Liu K, et al. Experimental study of a helium sorption cooler with low temperature fluctuation and long hold time below 1 K. Front Heat Mass Transf. 2025;23(3):739–50. doi:10.32604/fhmt.2025.064532. [Google Scholar] [CrossRef]

11. Duband L, Collaudin B. Sorption coolers development at CEA-SBT. Cryogenics. 1999;39(8):659–63. doi:10.1016/S0011-2275(99)00066-1. [Google Scholar] [CrossRef]

12. Klemencic GM, Ade PAR, Chase S, Sudiwala R, Woodcraft AL. A continuous dry 300 mK cooler for THz sensing applications. Rev Sci Instrum. 2016;87(4):045107. doi:10.1063/1.4945691. [Google Scholar] [PubMed] [CrossRef]

13. Kakac S, Bon B. A review of two-phase flow dynamic instabilities in tube boiling systems. Int J Heat Mass Transf. 2008;51(3–4):399–433. doi:10.1016/j.ijheatmasstransfer.2007.09.026. [Google Scholar] [CrossRef]

14. Marquardt ED, Le JP, Radebaugh R. Cryogenic material properties database. In: Ross RG, editor. Cryocoolers 11. Boston, MA, USA: Springer US; 2002. [Google Scholar]

15. Swartz ET, Pohl RO. Thermal boundary resistance. Rev Mod Phys. 1989;61(3):605–68. doi:10.1103/RevModPhys.61.605. [Google Scholar] [PubMed] [CrossRef]

16. Lei Y, Hong G, Quan J, Zhao YN, Li R, Wang G, et al. Design and development of a 4He sub-Kelvin sorption cooler. Cryogenics. 2024;141(2):103876. doi:10.1016/j.cryogenics.2024.103876. [Google Scholar] [CrossRef]

17. Chae DH, Won B, Choi J. Top-loading 4He sorption refrigerator. J Low Temp Phys. 2025;221(1):41–50. doi:10.1007/s10909-025-03313-3. [Google Scholar] [CrossRef]

18. Xi X, Wang J, Chen L, Zhou Y, Wang J. Progress and challenges of sub-kelvin sorption cooler and its prospects for space application. J Low Temp Phys. 2020;199(5):1363–81. doi:10.1007/s10909-020-02442-1. [Google Scholar] [CrossRef]

19. Pobell F. Matter and methods at low temperatures. 3rd ed. Berlin/Heidelberg, Germany: Springer; 2007. [Google Scholar]