This is a comment to John's answer, but since comments cannot have codeblocks, let me write a CW answer.

John's approach is to let $q_1,q_2,q_3,\dotsc$ be a list of proper prime powers so that $q_i - 1 \not\mid q_n-1$ for all $i < n$. Then, $$1 - \sum_i \frac{1}{(q_i-1)} + \sum_{i < j} \frac{1}{\mathrm{lcm}(q_i-1,q_j-1)} - \sum_{i < j < k} \frac{1}{\mathrm{lcm}(q_i-1,q_j-1,q_k-1)} \pm$$ approximates the natural density. Here is the SageMath code to calculate this sum:

def relevant_powers(limit):
    """lists the proper prime powers q such that q-1 is not a multiple of x-1 for previous powers"""
    powers = []
    for q in srange(2,limit+1):
        if q.is_prime_power() and not q.is_prime() and all([(q-1)%(x-1) > 0 for x in powers]):
            powers.append(q)
    return powers

def approximate_density(limit):
    """approximates the natural density of simple numbers using inclusion-exclusion principle"""
    sum = 0
    powers = relevant_powers(limit)
    for k in range(0,len(powers)+1):
        for selection in Subsets(powers,k):
            sum += (-1)^k * 1/lcm([x-1 for x in selection])
            print(float(sum)) # print progress
    return float(sum)

Then, executing approximate_density(1000000) yields 0.46211885256102014 at some point and it does not change anymore. (Of course, it does, but SageMath's precision is not enough.)

This is no proof, but this strongly indicates that there is a limit.

This is a comment to John's answer, but since comments cannot have codeblocks, let me write a CW answer.

John's approach is to let $q_1,q_2,q_3,\dotsc$ be a list of proper prime powers so that $q_i - 1 \not\mid q_n-1$ for all $i < n$. Then, $$1 - \sum_i \frac{1}{(q_i-1)} + \sum_{i < j} \frac{1}{\mathrm{lcm}(q_i-1,q_j-1)} - \sum_{i < j < k} \frac{1}{\mathrm{lcm}(q_i-1,q_j-1,q_k-1)} \pm$$ approximates the natural density. Here is the SageMath code to calculate this sum:

def relevant_powers(limit):
    """lists the proper prime powers q such that q-1 is not a multiple of x-1 for previous powers"""
    powers = []
    for q in srange(2,limit+1):
        if q.is_prime_power() and not q.is_prime() and all([(q-1)%(x-1) > 0 for x in powers]):
            powers.append(q)
    return powers

def approximate_density(limit):
    """approximates the natural density of simple numbers using inclusion-exclusion principle"""
    sum = 0
    powers = relevant_powers(limit)
    for k in range(0,len(powers)+1):
        for selection in Subsets(powers,k):
            sum += (-1)^k * 1/lcm([x-1 for x in selection])
            print(float(sum)) # print progress
    return float(sum)

Then, executing approximate_density(1000000) yields 0.46211885256102014 at some point and it does not change anymore. (Of course, it does, but SageMath's precision is not enough.)

However, executing approximate_density(5000000) yields 0.462118293064357 at some point and it does not change anymore. I am not sure what this means.