This paper introduces an efficient Monte Carlo (MC) algorithm for simulating "hardly-relaxing" systems, where multiple replicas with different temperatures are simultaneously simulated, and a virtual process exchanges configurations between replicas. This exchange process allows the system at low temperatures to escape from local minima. The algorithm is applied to the three-dimensional ±J Ising spin glass (SG) model, showing that it performs well even in the SG phase. The ergodicity time in this method is found to be much smaller than that of the multicanonical method. The time correlation function follows an exponential decay, with a relaxation time comparable to the ergodicity time at low temperatures, suggesting rapid relaxation through the exchange process. The method involves a compound system of many replicas, each with its own temperature. The system is simulated using a standard MC approach, with configurations of replicas exchanged based on energy cost. This allows replicas at low temperatures to escape local minima efficiently. The method ensures detailed balance, guaranteeing each replica equilibrates to its own canonical distribution. An advantage over simulated tempering is that no weighting factor is needed to equidistribute temperature probabilities. The algorithm is applied to the ±J Ising SG model, and the relaxation behavior is characterized by observing the ergodicity time and modified autocorrelation function. The autocorrelation function follows exponential decay, with a moderate relaxation time at low temperatures. The order parameter distribution P(q) is an even function, indicating equilibrium is reached. The method is organized as follows: Section 2 describes the exchange MC algorithm in detail. Section 3 discusses temperature determination. Section 4 reports the relaxational behavior of the method. Section 5 discusses the model and simulations, including temperature setting, relaxation, and order parameter distribution. The paper concludes with a discussion and summary, highlighting the efficiency of the method and its application to various models. The method is effective for systems with complex behavior, but may not work well for systems with first-order phase transitions. The paper is supported by numerical calculations and references to previous studies.This paper introduces an efficient Monte Carlo (MC) algorithm for simulating "hardly-relaxing" systems, where multiple replicas with different temperatures are simultaneously simulated, and a virtual process exchanges configurations between replicas. This exchange process allows the system at low temperatures to escape from local minima. The algorithm is applied to the three-dimensional ±J Ising spin glass (SG) model, showing that it performs well even in the SG phase. The ergodicity time in this method is found to be much smaller than that of the multicanonical method. The time correlation function follows an exponential decay, with a relaxation time comparable to the ergodicity time at low temperatures, suggesting rapid relaxation through the exchange process. The method involves a compound system of many replicas, each with its own temperature. The system is simulated using a standard MC approach, with configurations of replicas exchanged based on energy cost. This allows replicas at low temperatures to escape local minima efficiently. The method ensures detailed balance, guaranteeing each replica equilibrates to its own canonical distribution. An advantage over simulated tempering is that no weighting factor is needed to equidistribute temperature probabilities. The algorithm is applied to the ±J Ising SG model, and the relaxation behavior is characterized by observing the ergodicity time and modified autocorrelation function. The autocorrelation function follows exponential decay, with a moderate relaxation time at low temperatures. The order parameter distribution P(q) is an even function, indicating equilibrium is reached. The method is organized as follows: Section 2 describes the exchange MC algorithm in detail. Section 3 discusses temperature determination. Section 4 reports the relaxational behavior of the method. Section 5 discusses the model and simulations, including temperature setting, relaxation, and order parameter distribution. The paper concludes with a discussion and summary, highlighting the efficiency of the method and its application to various models. The method is effective for systems with complex behavior, but may not work well for systems with first-order phase transitions. The paper is supported by numerical calculations and references to previous studies.